Computer graphic system and method of multi-scale simulation of smoke

ABSTRACT

A multi-scale method is provided for computer graphic simulation of incompressible gases in three-dimensions with resolution variation suitable for perspective cameras and regions of importance. The dynamics is derived from the vorticity equation. Lagrangian particles are created, modified and deleted in a manner that handles advection with buoyancy and viscosity. Boundaries and deformable object collisions are modeled with the source and doublet panel method. The acceleration structure is based on the fast multipole method (FMM), but with a varying size to account for non-uniform sampling.

CROSS-REFERENCES TO RELATED APPLICATIONS

The present application is a non-provisional application of and claimsthe benefit and priority under 35 U.S.C. 119(e) of U.S. ProvisionalApplication No. 62/440,067, filed Dec. 29, 2016 entitled “MULTI-SCALEVORTICLE FLUID IMPROVEMENTS,” the entire content of which isincorporated herein by reference for all purposes.

BACKGROUND

Modeling the visual detail of a moving gas can be a laborious task. Theability to control sampling can vary from method to method. Foster andMetaxas [1997] solve the Navier-Stokes equation in a uniformly voxelizedgrid and further propose a popular unconditionally stable model. De Wittet al. [2012] represent the velocity using the Laplacian eigenvectors asa basis for incompressible flow. The Navier-Stokes equation can also besolved in a Lagrangian frame of reference with Smoothed ParticleHydrodynamics (SPH) [Gingold and Monaghan 1977], or with a VortexMethod, obtained by taking the curl of the Navier-Stokes equation[Cottet and Koumoutsakos 2000]. Vortex Methods can store data on points[Park and Kim 2005], curves [Angelidis et al. 2006b; Weissmann andPinkall 2010], or surfaces [Brochu et al. 2012], and can be solved withan integral, with a grid [Couet et al. 1981], or both [Zhang and Bridson2014]. Some original approaches don't use the Navier-Stokes equation[Elcott et al. 2007], but use Kelvin's theorem. Chern et al. [2016]solve the Schrodinger equation in a grid.

To model different characteristics at different resolutions, multipleapproaches can be combined and applied to a selective range of scales:Selle et al. [2005] advect vorticity carried by points and advectvorticity carried by sheets. Kim et al. [2008] advect procedurallygenerated detail. Pfaff et al. [2009] and Kim et al. [2012] show thatthe region of changing scale is located near boundaries, varying densityand temperature.

In the methods that compute pressure explicitly, the boundary conditioncan be met as part of computing the pressure. For the boundary conditionof Vortex Methods, a harmonic field may be needed, and the source anddoublet panel method can be the method of choice. Panel methods are wellresearched in aerodynamics. An introduction of panel methods is providedby Erickson [1990]. Both the source and the doublet term may berequired. Using only the source term may not prevent the flow frompassing through cup-shaped colliders, as opposed to spheroid collidersfor which the doublet contribution is very small. This may be the caseof the single layer potential of Zhang and Bridson [2014]. Using onlythe doublet term may not account for colliders with changing volume,since the doublet term is not capable of generating or consuming volume.Park and Kim [2005] use vorticle panels, which also do not handlecup-shaped colliders. Also, vorticle panels may not generate or consumevolume by themselves, and may induce a rotational field.

Therefore, there is a need for improved methods of computer graphicsimulation of an incompressible gas.

SUMMARY

Embodiments of the present invention provide a multi-scale method forcomputer graphic simulation of incompressible gases in three-dimensionswith resolution variation suitable for perspective cameras and regionsof importance. The dynamics may be derived from the vorticity equation.Lagrangian particles may be created, modified and deleted in a mannerthat handles advection with buoyancy and viscosity. Boundaries anddeformable object collisions may be modeled with the source and doubletpanel method. The acceleration structure may be based on the fastmultipole method (FMM), but with a varying size to account fornon-uniform sampling. Because the dynamics of the method is voxel free,one can freely specify the voxel resolution of the output density andvelocity while keeping the main shapes and timing.

These and other embodiments of the invention are described in detailbelow. For example, other embodiments are directed to systems, devices,and computer readable media associated with methods described herein.

A better understanding of the nature and advantages of embodiments ofthe present invention may be gained with reference to the followingdetailed description and the accompanying drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIGS. 1A-1B show a computer graphic renderings of a smoke trail thatstretches from near the camera (FIG. 1A) to the horizon (FIG. 1B), usinga single simulation according to some embodiments of the presentinvention.

FIGS. 2A-2C illustrate the fields of a vorticle: stream field (FIG. 2A),velocity field (FIG. 2B), and vorticity field (FIG. 2C), according tosome embodiments of the present invention.

FIG. 3 illustrates schematically a group of vorticles that arerelatively well separated from a target point p.

FIGS. 4A-4C and 5 illustrate a fast multipole Method (FMM) of evaluatingvelocity field of a plurality of vorticles according to some embodimentsof the present invention.

FIG. 6A illustrates 20 near vorticles that may induce a vector field.

FIGS. 6B-6E show pathlines induced by the near vorticles illustrated inFIG. 6A during one time step, for increasing numbers of time substepsusing straight segment lines, from 1 time substep to 3 time substeps(FIGS. 6B-6D), and to 64 time substeps (FIG. 6E), according to someembodiments of the present invention.

FIGS. 6F-6H show pathlines induced by the near vorticles in FIG. 6A forincreasing numbers of time substeps, using twists constructed using theLie group SE(3), from 1 time substep to 3 time substeps (FIGS. 6F-6H),according to some embodiments of the present invention.

FIGS. 7A-7C illustrate the three main steps of a fast multipole method(FMM) that may be used to compute the near coefficients of a faradvected field's near expansion according to some embodiments of thepresent invention.

FIG. 8 illustrates schematically a plurality of vorticles of varioussizes that may be generated in a virtual space surrounding a cameralocated at p according to some embodiments of the present invention.

FIG. 9 illustrates an exemplary octree cell structure for storingvorticles according to some embodiments of the present invention.

FIGS. 10A-10B illustrate how individual vorticles inside the box (FIG.10A) are transformed into far coefficients relative to the box's center(FIG. 10B) for evaluations outside of the 26 neighboring cells shownwith a dashed line according to some embodiments of the presentinvention.

FIGS. 11A-11B illustrate how the far field of a cell (FIG. 11A) istranslated to the parent cell (FIG. 11B) according to some embodimentsof the present invention.

FIGS. 12A-12B illustrate how, to evaluate the far field of a cell insidethe dashed red line (FIG. 12A), the far field is transformed into a nearfield (FIG. 12B) according to some embodiments of the present invention.

FIGS. 13A-13B illustrate how the near field of a cell (FIG. 13A) istranslated to a child cell (FIG. 13B) according to some embodiments ofthe present invention.

FIGS. 14A-14B illustrate an initial position of a boundary object (e.g.,a tea pot) moving into points (the lattice) (FIG. 14A), and how theharmonic field may warp the space in an incompressible and irrotationalmanner with a slip boundary (FIG. 14B) according to some embodiments ofthe present invention.

FIGS. 15A-15C illustrate a first vorticle (FIG. 15A) is stretched into asecond vorticle (FIG. 15B), where the mean energy is preserved while theenstrophy changes, and the second vorticle is then resampled in anunstretched third vorticle (FIG. 15C) with the same mean energy andenstrophy as the second vorticle, according to some embodiments of thepresent invention.

FIG. 16 illustrates a density particle emits a ring of vorticles (red)according to some embodiments of the present invention.

FIG. 17A illustrates a constant flow moves at 10 ms⁻¹ with v=0.1 arounda static pole of radius 1 (red): the surface sheds vorticles (blackarrows), according to some embodiments of the present invention.

FIG. 17B illustrates a slice of the vorticity induced by the vorticlesthat reveals the emergent behavior of a von Karman vortex streetaccording to some embodiments of the present invention.

FIGS. 18A-18E illustrate schematically a method of forming density voxelgrids according to some embodiments of the present invention.

FIG. 19 is a flowchart illustrating a method of computer graphic imageprocessing according to some embodiments of the present invention.

FIGS. 20A-20B illustrate a simulated smoke column according to someembodiments of the present invention.

FIGS. 21A-21D illustrate an example of a simulator that may takeadvantage of a camera frustum according to some embodiments of thepresent invention.

FIG. 22A illustrates how tiny vorticles (shown in white) may interactwith large vorticles (shown in black) according to some embodiments ofthe present invention.

FIG. 22B shows the density advected using the vorticles illustrated inFIG. 22A.

FIGS. 23A-23B illustrate how the harmonic field may play the same roleas the pressure term in the classic grid-based Navier-Stokes equationaccording to some embodiments of the present invention.

FIG. 24 is a simplified block diagram of a system for creating computergraphics imagery (CGI) and computer-aided animation that may implementor incorporate various embodiments of the present invention.

FIG. 25 is a block diagram of a computer system or informationprocessing device that may incorporate an embodiment, be incorporatedinto an embodiment, or be used to practice any of the innovations,embodiments, and/or examples found within this disclosure.

FIG. 26A shows two almost exactly overlapping curves of the largestcomponent of the left (red) and right (blue) side of Eq. (54) accordingto some embodiments of the present invention; FIG. 26B shows the densityfield ρ_(j) for r_(j)=1 and m_(j)=1 according to some embodiments of thepresent invention.

FIGS. 27A-27D illustrate that the integral over any triangle can bedecomposed into 6 right-angle triangle integrals according to someembodiments of the present invention.

DETAILED DESCRIPTION

Embodiments of the present invention provide a multi-scale method forcomputer graphic simulation of incompressible gases in three-dimensionswith resolution variation suitable for perspective cameras and regionsof importance. The dynamics may be derived from the vorticity equation.Lagrangian particles may be created, modified and deleted in a mannerthat handles advection with buoyancy and viscosity. Boundaries anddeformable object collisions may be modeled with the source and doubletpanel method. The acceleration structure may be based on the fastmultipole method (FMM), but with a varying size to account fornon-uniform sampling. Because the dynamics of the method is voxel free,one can freely specify the voxel resolution of the output density andvelocity while keeping the main shapes and timing.

The FMM was developed by Geengard and Rokhlin [1987] to solve the N-bodyproblem in linear complexity. An introduction of the FMM is provided byYing [2012]. Unlike hybrid methods, the methods disclosed herein canhandle all ranges of scales in one model. Unlike SPH where small andlarge scales cannot overlap, the samples in these methods can overlapfor all scales. Unlike grid-based methods that have a uniform samplingfor the finest resolution, these methods can handle a sparse structurefor all scales. Embodiments may provide: (1) a hierarchical and sparsemethod for free flow advection; (2) a hierarchical and sparse method forboundaries; (3) an unconditionally stable method for vorticitystretching; and (4) a new model for buoyancy.

FIGS. 1A-1B show a computer graphic renderings of a smoke trail thatstretches from near the camera (FIG. 1A) to the horizon (FIG. 1B), usinga single simulation according to some embodiments of the presentinvention.

A. Dynamic Model

To obtain the equation of dynamics, the curl operator may be applied tothe following form of the Navier-Stokes equation of a Newtonian fluid:

$\begin{matrix}{{\rho \frac{d\overset{\rightarrow}{u}}{dt}} = {{\mu {\nabla^{2}\overset{\rightarrow}{u}}} + {\rho \overset{\rightarrow}{F}} - {\nabla p}}} & (1)\end{matrix}$

where {right arrow over (u)} is velocity, t is time, ρ is density, μ isviscosity, {right arrow over (F)} is external force, and p is pressure.

Assuming that μ is constant, dividing both sides of Eq. (1) by ρ,replacing

$\frac{\mu}{\rho}$

with ν, and use the identity ∇ρ/ρ=∇ log(ρ), then the following equationmay be obtained:

$\begin{matrix}{{\frac{d\overset{\rightarrow}{\omega}}{dt} = {{\left( {\overset{\rightarrow}{\omega} \cdot \nabla} \right)\overset{\rightarrow}{u}} + {v\; {\nabla^{2}\overset{\rightarrow}{\omega}}} + {{\nabla{\log (\rho)}} \times \left( {\overset{\rightarrow}{F} - \frac{d\overset{\rightarrow}{u}}{dt}} \right)} + {\nabla{\times \overset{\rightarrow}{F}}}}},} & (2)\end{matrix}$

where the vorticity {right arrow over (ω)} is the curl of the velocity{right arrow over (u)}. This is the complete equation, it includesboundaries and all forces. Eq. (2) means that the vorticity {right arrowover (ω)} evolves over time by advecting particles carrying {right arrowover (ω)}, and by stretching {right arrow over (ω)} according to thevelocity {right arrow over (u)}, with kinematic viscosity ν, buoyancyand boundary interaction specified by density ρ and external force{right arrow over (F)}:

$\begin{matrix}{{\overset{\rightarrow}{F}(p)} = \left\{ {\begin{matrix}{\overset{\rightarrow}{g} + \overset{\rightarrow}{e}} & {{if}\mspace{14mu} p\mspace{14mu} {outside}\mspace{14mu} {solid}\mspace{14mu} {objects}} \\\overset{\rightarrow}{f} & {otherwise}\end{matrix},} \right.} & (3)\end{matrix}$

where {right arrow over (g)} is the constant for gravity, {right arrowover (e)} is a set of user defined external forces, and {right arrowover (f)} is the acceleration at the objects' boundaries, suitable fordeformable objects. The term “advecting” may refer to “moving” (e.g.,particles). Solving {right arrow over (ω)} with Eq. (2) does not involvethe pressure term. Instead, the velocity {right arrow over (u)} isobtained from {right arrow over (ω)} in two parts: an advected field{right arrow over (v)} derived from {right arrow over (ω)}, and aharmonic field {right arrow over (h)} derived from {right arrow over(v)} at the boundaries:

{right arrow over (u)}={right arrow over (v)}+{right arrow over(h)}  (4)

The fields {right arrow over (v)} and {right arrow over (h)} satisfy thefollowing identities:

∇·{right arrow over (v)}=0 ∇·{right arrow over (h)}=0 ∇×{right arrowover (h)}=0  (5)

B. Advected Field

The advected field {right arrow over (v)} is the purely rotational anddynamic portion of {right arrow over (u)} that can be computed directlyfrom {right arrow over (ω)}. It is defined with the Biot-Savart law:

$\begin{matrix}{{\overset{\rightarrow}{v}(p)} = {\frac{1}{4\pi}{\int_{x \in R^{3}}{\frac{{\overset{\rightarrow}{\omega}(x)} \times \left( {p - x} \right)}{{{p - x}}^{3}}{dx}}}}} & (6)\end{matrix}$

Eq. (6) may mean that {right arrow over (v)} is induced by a continuumof weighted rotations singular at p=x. To remove the singularity, adiscrete overlapping partitioning may be used, where each element i is avorticle. A vorticle may also be referred to as a vortex particle or avortex blob.

A vorticle may be defined by an integrable vorticity field {right arrowover (ω)}_(i) centered at p_(i):

$\begin{matrix}{{\overset{\rightarrow}{v}(p)} = {{\frac{1}{4\pi}{\int_{x \in R^{3}}{\frac{{\overset{\rightarrow}{\omega}(x)} \times \left( {p - x} \right)}{{{p - x}}^{3}}{dx}}}} = {\Sigma_{i}{{{\overset{\rightarrow}{v}}_{i}\left( {p - p_{i}} \right)}.}}}} & (7)\end{matrix}$

A motivation to use particles instead of a topologically connectedstructure, such as filaments, may be to avoid the rapid entanglementoccurring around objects and other filaments. A vorticle i may beassociated with a vorticity field {right arrow over (ω)}_(i), a velocityfield {right arrow over (v)}_(i) and a stream field {right arrow over(ψ)}_(i). The following explicit expressions may be chosen for thestream field {right arrow over (–)}_(i), the velocity field {right arrowover (v)}_(i), and the vorticity field {right arrow over (ω)}_(i), using

${{\xi_{i}(l)} = \frac{1}{\sqrt{r_{i}^{2} + l^{2}}}},$

size r_(i) and rotation strength {right arrow over (w)}_(i):

{right arrow over (ψ)}_(i)({right arrow over (x)})=ξ_(i)(∥{right arrowover (x)}∥){right arrow over (w)} _(i)

{right arrow over (v)}_(i)({right arrow over (x)})=ξ_(i)(∥{right arrowover (x)}∥)³ {right arrow over (w)} _(i) ×{right arrow over (x)}

{right arrow over (ω)}_(i)({right arrow over (x)})=ξ_(i)(∥{right arrowover (x)}∥)³(2{right arrow over (w)} _(i)−3ξ_(i)(∥{right arrow over(x)}∥)² {right arrow over (x)}×({right arrow over (w)} _(i) ×{rightarrow over (x)}))  (8)

FIGS. 2A-2C illustrate the fields of a vorticle: stream field (FIG. 2A),velocity field (FIG. 2B), and vorticity field (FIG. 2C) according tosome embodiments, where {right arrow over (w)}_(i) is up. The fieldsshown in FIGS. 2A-2C satisfy the following identities:

∇×{right arrow over (ψ)}_(i)={right arrow over (ν)}_(i) ∇×{right arrowover (v)}_(i)={right arrow over (ω)}_(i) ∇·{right arrow over (ν)}_(i)=0.

These functions are chosen because they are non-singular since r_(i)>0,they are smooth across {right arrow over (x)}=0, and because thevorticle's rotation strength {right arrow over (w)}_(i) is asymptotic toits vorticity {right arrow over (ω)}_(i) when decreasing r_(i). Also, avorticle's mean energy E_(i) can have a closed form:

$\begin{matrix}{{E_{i} = {{\frac{1}{2}{\int{{{{\overset{\rightarrow}{v}}_{i}\left( \overset{\rightarrow}{x} \right)}}^{2}d\overset{\rightarrow}{x}}}} = \frac{\pi^{2}{{\overset{\rightarrow}{w}}_{i}}^{2}}{4\; r_{i}}}},} & (9)\end{matrix}$

which may be convenient for splitting a vorticle emitter's energy Eamong N emitted vorticles:

$\begin{matrix}{{\overset{\rightarrow}{w}}_{i} = {\frac{2}{\pi}\sqrt{\frac{{Er}_{i}}{N}}{X.}}} & (10)\end{matrix}$

The variable X is a random vector on the unit sphere, or a normalizedcurl of an input field. A set of vorticles induces a stream field {rightarrow over (ψ)}, an advected field {right arrow over (v)} and avorticity field {right arrow over (ω)}:

{right arrow over (ψ)}(p)=Σ{right arrow over (ψ)}_(i)(p−p _(i))

{right arrow over (v)}(p)=Σ{right arrow over (v)} _(i)(p−p _(i))

{right arrow over (ω)}(p)=Σ{right arrow over (ω)}_(i)(p−p _(i))  (11)

Therefore the following identities hold true by construction:

∇×{right arrow over (ψ)}={right arrow over (v)} ∇×{right arrow over(v)}={right arrow over (ω)} ∇·{right arrow over (v)}=0 ∇·{right arrowover (ω)}=0

Multi-Scale Fast Multipole Method (FMM)

The cost of evaluating {right arrow over (v)} at every vorticle centerp_(i) directly with Eq. (11) scales quadratically with the number ofvorticles, and may be prohibitive for large numbers of vorticles. Thefast multipole method (FMM) may provide a way to evaluate {right arrowover (v)}(p) using only the near cells of p in an octree. The FMM is anumerical technique that was developed to speed up the calculation oflong-ranged forces in a N-body problem. It accomplishes this by using amultipole expansion, which allows one to group sources that lie closetogether and treat them as if they are a single source.

FIG. 3 illustrates schematically a group of vorticles 310 that arerelatively well separated from a target point p. According to an FMMmethod in some embodiments, when sufficiently far away from p, the groupof vorticles 310 may be approximated as a single source centered on acell 320.

FIGS. 4A-4C and 5 illustrate schematically an FMM method of evaluatingvelocity field of a plurality of vorticles according to someembodiments. Assume that vorticles 410 are approximately uniformlydistributed in a three-dimensional space, as illustrated in FIG. 4A. Thethree-dimensional space may be divided into cubical cells 420. For atarget point located sufficiently far away from a cell 420, thevorticles 410 inside the cell 420 may be approximated as a single source430 positioned at the center of the cell 420, as illustrated in FIG. 4B.

An octree cell structure may be constructed to include a plurality ofoctree levels. Traversing the octree up to a parent level, eight nearestneighbor child cells 420 (e.g., four in the plane of the paper and fourin a plane below the plane of the paper) may be approximated as a singlesource 450 positioned at the center of the parent cell 440, asillustrated in FIG. 4C. This process can continue to the next parentlevel in a similar fashion.

In some embodiments, a cell that contains the camera position p and thecells that are direct neighbors of that cell may be referred to as nearcells. All other cells may be referred to as far cells. For instance, inthe example illustrated in FIG. 5, the cells included in thethick-outlined box 510 surrounding the camera position p may be referredto as near cells, and the cells outside the thick-outlined box 510 maybe referred to as far cells. In some embodiments, the vorticles in thecell and in the 1-cell ring of p may be referred to as near vorticles,and all other vorticles may be referred to as far vorticles. Note thatthe octree may store vorticles of different sizes, and there arepotentially near and far vorticles at multiple levels.

The advected field may be split into near and far fields:

{right arrow over (v)}={right arrow over (v)} _(near) +{right arrow over(v)} _(far).  (12)

1. Near Advected Field

The near advected field may be evaluated by summing the velocity fieldsof the near vorticles, in the cell and 1-cell ring across all levels:

$\begin{matrix}{{{\overset{\rightarrow}{v}}_{near}(p)} = {\sum\limits_{i \in {{near}{(p)}}}{{\overset{\rightarrow}{v}}_{i}\left( {p - p_{i}} \right)}}} & (13)\end{matrix}$

The near field, however, can have a very large mean energy, especiallynear user defined vorticle emitters. Also, if r_(i) approaches 0, then{right arrow over (v)}_(i)(p_(i)) becomes infinite. For the simulationto be tolerant of any input, an advection method that handles the mostextremely coiled trajectory in a single time step may be used. For that,the velocity may be replaced with a displacement constructed from theunion of twists, similar to the method of Angelidis et al. [2006a]. Atwist is an element of the Lie group SE(3). A twist may be constructedby combining the rotations induced by vorticles. Lie group of twistswith vorticles may be used because it is known that the motion isderived from rotations. The rotations may be expressed in the Liealgebra se(3), and summed. An exponential map may be used to obtain thetwist in SE(3). The following may be obtained:

$\begin{matrix}{\mspace{79mu} {{{\overset{\rightarrow}{k}}_{0} = {\Delta_{t}{\sum\limits_{i \in {{near}{(p)}}}{{{\overset{\rightarrow}{\xi}}_{i}\left( {p - p_{i}} \right)}^{2}{{\overset{\rightarrow}{\psi}}_{i}\left( {p - p_{i}} \right)}}}}}\mspace{20mu} {{\overset{\rightarrow}{k}}_{1} = {\Delta_{t}{\sum\limits_{i \in {{near}{(p)}}}{{\overset{\rightarrow}{v}}_{i}\left( {p - p_{i}} \right)}}}}{{{\overset{\rightarrow}{v}}_{near}(p)} = {\frac{1}{\Delta_{t}}\left( {{\overset{\rightarrow}{k}}_{1} + {{\overset{\rightarrow}{k}}_{0} \times \left( {{\frac{{{\overset{\rightarrow}{k}}_{0}} - {\sin \left( {{\overset{\rightarrow}{k}}_{0}} \right)}}{{{\overset{\rightarrow}{k}}_{0}}^{3}}{\overset{\rightarrow}{k}}_{0} \times {\overset{\rightarrow}{k}}_{1}} + {\frac{1 - {\cos \left( {{\overset{\rightarrow}{k}}_{0}} \right)}}{{{\overset{\rightarrow}{k}}_{0}}^{2}}{\overset{\rightarrow}{k}}_{1}}} \right)}} \right)}}}} & (14)\end{matrix}$

This method may underestimate the length of the trajectory, butconverges when decreasing the time step Δ_(i) and guarantees stability.It is observed that strong vorticles induce extremely coiled pathlinesin a single time step similar to the pathlines integrated with tinysteps, and also that a small number of vorticles placed on a ring movesalong a straight and stable trajectory.

FIG. 6A illustrates 20 near vorticles that may induce a vector field.FIGS. 6B-6E show pathlines induced by the near vorticles illustrated inFIG. 6A during one time step, for increasing numbers of time substepsusing straight segment lines (referred to as linear time steps), from 1time substep to 3 time substeps (FIGS. 6B-6D), and to 64 time substeps(FIG. 6E), according to some embodiments of the present invention.Points are scattered with random colors. As illustrated, the quality ofthe advection increases rather slowly when taking linear time steps.FIGS. 6F-6H show pathlines induced by the near vorticles in FIG. 6A forincreasing numbers of time substeps, using twists constructed using theLie group SE(3), from 1 time substep to 3 time substeps (FIGS. 6F-6H),according to some embodiments of the present invention. As illustrated,using the Lie group of twists, thousands of rotations can be done in asingle time step.

2. Far Advected Field

The far advected field's near expansion may be defined using the nearcoefficients {right arrow over (κ)}_(lmn) of the deepest cell containingp:

{right arrow over (v)} _(far)(p)=Σ_(l=0) ^(l) ^(max) Σ_(m=−l)^(l)Σ_(n=1) ^(n) ^(max) ∇(Y _(lm)(p)∥p∥ ^(n))×{right arrow over(κ)}_(lmn)  (15)

The function Y_(lm)(p) denotes the real spherical harmonics:

$\begin{matrix}{{K_{lm} = \sqrt{\frac{\left( {{2l} + 1} \right){\left( {l - {m}} \right)!}}{4{{\pi \left( {l + {m}} \right)}!}}}}{{Y_{lm}(p)} = {K_{lm}{P_{l,{m}}\left( {\cos (\varphi)} \right)}\left\{ \begin{matrix}{\sqrt{2}{\cos \left( {m\; \theta} \right)}} & {{{if}\mspace{14mu} m} > 0} \\1 & {{{if}\mspace{14mu} m} = 0} \\{\sqrt{2}{\sin \left( {{- m}\; \theta} \right)}} & {{{if}\mspace{14mu} m} < 0}\end{matrix} \right.}}} & (16)\end{matrix}$

The function P_(lm) denotes the associated Legendre polynomials, withCartesian coordinates mapping to spherical coordinates as

$\left( {\theta,\varphi,r} \right) = {\left( {{a\; \tan \; 2\left( {p_{y},p_{x}} \right)},{a\; {\cos\left( \frac{p_{z}}{p} \right)}},{p}} \right).}$

The first few basis functions are given in Appendix 2.

FIGS. 7A-7C illustrate the three main steps of the FMM that may be usedto compute the near coefficients {right arrow over (κ)}_(lmn). First, anoctree can be built where the near coefficients {right arrow over(κ)}_(lmn) of all cells are initialized to 0. Level 0 denotes theoctree's root. The 189 cells that are children of the parent's neighborsbut not direct neighbors may be referred to as far interacting cells. Acell's far coefficient may be computed by translating the 8 children'sfar coefficient (FIG. 7A). A cell's near coefficient may be computed bytranslating the parent's near coefficient (FIG. 7B), and transformingthe 189 far interacting cells' far coefficient (FIG. 7C).

To make the algorithm multiscale, vorticles of various sizes r_(i) maybe generated in a three-dimensional space according to some embodiments,as illustrated in FIG. 8. For examples, the vorticles may include smallvorticles 810 and 812, medium vorticles 820, medium large vorticles 830,and large vorticles 840, 842, and 844. In FIG. 8, the point p denotesthe camera position. In some embodiments, the size of each vorticle mayincrease as the distance from the vorticle to the camera position pincreases. In some other embodiments, vorticles of various sizes mayoccupy a same spatial region. For instance, in the example illustratedin FIG. 8, the small vorticle 812 may be located next to the largevorticles 842 and 844.

The vorticles may be stored in the internal cells of size

$a = {2^{\lceil\frac{\log \mspace{11mu} r_{i}}{{\log \mspace{11mu} 2}\;}\rceil}.}$

FIG. 9 illustrates an exemplary cell structure. In some embodiments, thecell size a may be expressed as a=2^(−L), where

$L = \frac{\log \mspace{11mu} r_{i}}{{\log \mspace{11mu} 2}\;}$

is the level in an octree. When traversing up the octree from a childcell to a parent cell (i.e., from lower level L to higher level L−1),larger vorticles may be collected.

According to some embodiments, the steps in an FMNI method may include:

(1) Compute the far coefficients {right arrow over (μ)}_(lmn) of allcells with vorticles, as described further below and illustrated inFIGS. 10A-10B. The far interacting cells may also be allocated all theway to the root, so that steps (2) and (3) may propagate validcoefficients to the deepest cell, to cover all positions with theoctree.

(2) Traverse the tree from the leaves to the root level 0, translatingand accumulating the children's far coefficients {right arrow over(μ)}_(lmn) into the parent's far coefficient {right arrow over(μ)}_(l′m′n′)′, as described further below and illustrated in FIGS.11A-11B.

(3) Traverse the tree from the root level 0 to the leaves: first, theparent's near coefficient is translated, as described further below andillustrated in FIGS. 12A-12B; second, the far interacting cells' farcoefficient is transformed, as described further below and illustratedin FIGS. 13A-137B. Each far interacting cell's set of coefficients{right arrow over (κ)}_(lmn). is added to the set of coefficients storedin the current cell.

(4) After step (3), the coefficients at the leaves represent all farvorticles are then used in Eq. (15).

a) Vorticle Far Coefficients

The far advected field may be decomposed around the origin using the farcoefficients {right arrow over (μ)}_(lmn):

$\begin{matrix}{{{\overset{\rightarrow}{v}}_{far}(p)} = {\sum\limits_{l = 0}^{l_{\max}}{\sum\limits_{m = {- l}}^{l}{\sum\limits_{n = 1}^{n_{\max}}{{\nabla\left( {{Y_{lm}(p)}\frac{1}{{p}^{n}}} \right)} \times {\overset{\rightarrow}{\mu}}_{lmn}}}}}} & (17)\end{matrix}$

Unlike a traditional expansion, an advantage of this expansion may bethat it includes the vorticles' size. To compute the far coefficients{right arrow over (μ)}_(lmn), a far expansion of a single vorticle'sstream function may be made as:

$\begin{matrix}{{{\overset{\rightarrow}{\psi}}_{i}\left( {p - p_{i}} \right)} = {{\overset{\rightarrow}{w}}_{i}{\xi_{i}\left( {p} \right)}{\sum\limits_{k}{{P_{k}\left( {\frac{p \cdot p_{i}}{{p}{p_{i}}},{\tau_{i}\left( {p} \right)}} \right)}\frac{{p_{i}}^{k}}{{p}^{k}}}}}} & (18)\end{matrix}$

This can be achieved by identifying the terms of a Taylor expansion nearϵ=0 of ξ_(i)(∥p−p_(i)∥)=ξ_(i)(∥p∥)(1+τ_(i)(∥p∥ϵ)^(−1/2), where

${{\tau_{i}(l)} = \frac{1}{1 + {r_{i}^{2}/l^{2}}}},$

and where P_(k)(z, τ) is a modified Legendre polynomial, that is, aLegendre polynomial where the ith non-zero coefficient of the largestpower of z is multiplied by τ_(i) ^(k−i). The first few are:

$\begin{matrix}{{{P_{0}\left( {\alpha,\tau} \right)} = 1}{{P_{1}\left( {\alpha,\tau} \right)} = {\alpha\tau}}{{P_{2}\left( {\alpha,\tau} \right)} = {{- \frac{\tau}{2}} + \frac{3\alpha^{2}\tau^{2}}{2}}}} & (19)\end{matrix}$

This modified Legendre decomposition may be projected on the Legendrepolynomials:

$\begin{matrix}{{{\overset{\rightarrow}{\psi}}_{i}\left( {p - p_{i}} \right)} = {{\overset{\rightarrow}{w}}_{i}{\sum\limits_{l}{{P_{l}\left( \frac{p \cdot p_{i}}{{p}{p_{i}}} \right)}{B_{l}\left( {{p},{p_{i}}} \right)}\mspace{14mu} {where}\text{:}}}}} & (20) \\{{{B_{l}\left( {r,d_{i}} \right)} = {{\xi_{i}(r)}{\sum\limits_{k = l}^{\infty}{{\beta_{kl}\left( {\tau_{i}(r)} \right)}\frac{d_{i}^{k}}{r^{k}}}}}}{{\beta_{kl}(\tau)} = {\int_{- 1}^{1}{\frac{{2l} + 1}{2}{P_{l}(\alpha)}{P_{k}\left( {\alpha,\tau} \right)}d\; \alpha}}}} & (21)\end{matrix}$

A series expansion in 1/r may then be perform to obtain the radialcoefficients b_(ln):

$\begin{matrix}{B_{l} = {\sum\limits_{n = 1}^{\infty}{b_{\ln}\frac{1}{r^{n}}}}} & (22)\end{matrix}$

Finally, the Legendre polynomial may be projected on the sphericalharmonics to obtain the spherical coefficients y_(lm):

$\begin{matrix}{y_{lm} = {\int_{\theta = 0}^{2\pi}{\int_{\varphi = 0}^{\pi}{{P_{l}\left( \frac{p_{r\; {\theta\varphi}} \cdot p_{i}}{p_{i}} \right)}{Y_{lm}\left( {\theta,\varphi} \right)}d\; \theta \; d\; \varphi}}}} & (23)\end{matrix}$

Thus the far coefficients {right arrow over (μ)}_(lmn) of one vorticleis assembled from the radial coefficient, the spherical coefficient, andthe vorticle's rotation strength {right arrow over (w)}_(i):

{right arrow over (μ)}_(lmn)=y_(lm)b_(ln){right arrow over(w)}_(i)  (24)

The first few far coefficients are given in Appendix 3. The far field ofmultiple vorticles may be obtained simply by summing the coefficients ofeach vorticle.

FIGS. 10A-10B illustrates how individual vorticles inside the box (FIG.10A) are transformed into far coefficients relative to the box's center(FIG. 10B) for evaluations outside of the 26 neighboring cells shownwith a dashed line. This reduction of complexity to a basis of constantsize makes this method scalable. For evaluating {right arrow over (v)}outside of the octree, the far expansion of Eq. (17) is used at theroot.

b) Far Coefficients Translation

FIGS. 11A-11B illustrates how the far field of a cell (FIG. 11A) istranslated to the parent cell (FIG. 11B). The matching area is outlinedin red. The 26 neighboring cells are shown with a dashed line. Tocompute the parent's far coefficients {right arrow over (μ)}_(l′m′n′)from child coefficients {right arrow over (μ)}_(lmn) relative to offsetlocation c, a series expansion of {right arrow over (ψ)}_(far) (p+c) in1/r can be made and projected on the translated basis Y_(l′m′). Thefirst few are given in Appendix 4.

Since the relative spatial configuration between parent and childrennodes are repeated in the octree, the translation coefficients can beprecomputed, and the translation becomes a mere sum of weighted vectors.This is also the case for all other coefficients.

c) Far-to-Near Coefficients Transform

FIGS. 12A-12B illustrate how, to evaluate the far field of a cell insidethe dashed red line (FIG. 12A), the far field is transformed into a nearfield (FIG. 12B). For far interacting cells, the far expansioncoefficients {right arrow over (μ)}_(lmn) are transformed into the nearexpansion coefficients {right arrow over (κ)}_(l′m′n′). To compute thenear coefficients {right arrow over (κ)}_(l′m′n′) of the translated farcoefficient {right arrow over (μ)}_(lmn) to a new location c, a seriesexpansion of {right arrow over (ψ)}_(far)(p+c) in r may be made andprojected on the translated basis Y_(l′m′) to obtain the coefficient ofthe translated field. Since some of the far coefficients are zero, theysimplify, and the first few are given in Appendix 5.

d) Near Coefficients Translation

The coefficients from parent nodes {right arrow over (κ)}_(lmn) can bepropagated to the coefficient of children nodes {right arrow over(κ)}′_(l′m′n′), by translating the parent's coefficient to each child'scenter, and adding the coefficient. To compute the coefficients {rightarrow over (κ)}′_(l′m′n′) of the translated coefficient {right arrowover (κ)}_(lmn) to a new location c, a series expansion of {right arrowover (ψ)}_(near)(p+c) in r can be made and projected on the translatedbasis Y_(l′m′) to obtain the coefficient of the translated field. Theysimplify and the first few are given in Appendix 6.

FIGS. 13A-13B illustrate how the near field of a cell (FIG. 13A) istranslated to a child cell (FIG. 13B).

C. Harmonic Field

The harmonic field {right arrow over (h)} in Eq. (4) can be computedusing the source and doublet panel method. In contrast with the couplingof grid based solvers, which solve for both pressure and boundariessimultaneously in a volume domain, the source and doublet panel methodsolves the boundaries' degrees of freedom on a surface domain, with asolution smoothly defined between the boundaries and infinity. Given amanifold surface defined by n triangle panels δΩ_(j) with normal {rightarrow over (n)}_(j) pointing outside, field {right arrow over (h)} isdefined as:

{right arrow over (h)}=Σ _(j)σ_(j) ∇S _(j)+μ_(j) ∇D _(j).  (25)

Each triangle induces a source panel field S_(j) and a doublet panelfield D_(j), defined by surface integrals:

$\begin{matrix}{{{S_{j}(p)} = {\int_{\;}^{\;}{\int_{{\delta\Omega}_{j}}^{\;}{\frac{1}{{p - x}}{dx}}}}}{{D_{j}(p)} = {- {\int_{\;}^{\;}{\int_{{\delta\Omega}_{j}}^{\;}{{{\overset{\rightarrow}{n}}_{j} \cdot \frac{p - x}{{{p - x}}^{3}}}{dx}}}}}}} & (26)\end{matrix}$

The gradients of the source and the doublet are both irrotational andincompressible, and therefore ∇·{right arrow over (h)}=0 and ∇×{rightarrow over (h)}=0. To find σ_(j) and μ_(j) such that field {right arrowover (h)} cancels field {right arrow over (v)} across all boundaries,the linear system of the direct approach can be solved by placing acontrol point p_(i) at the center of each panel. This produces a linearsystem of 2n equations with 2n variables:

$\begin{matrix}\left\{ \begin{matrix}{0 = {{\sum\limits_{j}{\sigma_{j}{S_{j}\left( p_{i} \right)}}} + {\mu_{j}{D_{j}\left( p_{i} \right)}}}} \\{{{\overset{\rightarrow}{n}}_{i} \cdot \left( {{\overset{\rightarrow}{s}\left( p_{i} \right)} - {\overset{\rightarrow}{v}\left( p_{i} \right)}} \right)} = {{\sum\limits_{j}{\sigma_{j}{{\overset{\rightarrow}{n}}_{i} \cdot {\nabla{S_{j}\left( p_{i} \right)}}}}} + {\mu_{j}{{\overset{\rightarrow}{n}}_{i} \cdot {\nabla{D_{j}\left( p_{i} \right)}}}}}}\end{matrix} \right. & (27)\end{matrix}$

where {right arrow over (s)}(p_(i)) the velocity at the center of paneli. If the boundary isn't moving, then {right arrow over (s)}(p_(i))=0.For numerical evaluations, analytical integrals for S_(j) and ∇S_(j),D_(j) and ∇D_(j) may be used, provided in Appendices 9 and 10.

For evaluating the harmonic field {right arrow over (h)}, the harmonicfield in near and far fields may be split:

{right arrow over (h)}={right arrow over (h)} _(near) +{right arrow over(h)} _(far).  ((28)

The panels may be stored in an octree structure similar to the one usedfor the advected field. Note that this can be a second octree, sincevorticles sample space between colliders, whereas the panels sample thesurface of the colliders, and there is generally no cell correlationbetween the two.

1. Near Harmonic Field

The near harmonic field may be evaluated directly by summing theharmonic fields of the near panels, in the cell and 1-cell ring acrossall levels, in a manner similar to Eq. (13):

$\begin{matrix}{{{\overset{\rightarrow}{h}}_{near}(p)} = {{\sum\limits_{j \in {{near}{(p)}}}{\sigma_{j}{\nabla{S_{j}(p)}}}} + {\mu_{j}{\nabla{D_{j}(p)}}}}} & (29)\end{matrix}$

FIGS. 14A-14B illustrate an initial position of a boundary object (e.g.,a tea pot) moving into points (the lattice) (FIG. 14A), and how theharmonic field may warp the space in an incompressible and irrotationalmanner with a slip boundary (FIG. 14B).

2. Far Harmonic Field

The far harmonic field is similar to the far advected field, but insteadof the curl of a field with vector coefficients, it is the gradient of afield with scalar coefficients. The near harmonic coefficients κ_(lmn)can be computed using the same FMNI algorithm, but only with1-dimensional coefficients, and with samples on the surface:

{right arrow over (h)} _(far)(p)=Σ_(l=0) ^(l) ^(max) Σ_(m=−l)^(l)Σ_(n=1) ^(n) ^(max) κ_(lmn)∇(Y _(lm)(p)∥p∥ ^(n))  (30)

Eq. (30) is similar to Eq. (15). To derive the far coefficients, eachsource and doublet panel of area a_(i) may be approximated with a point:

$\begin{matrix}{{\overset{\rightarrow}{h}}_{far} = {\nabla\left( {{{\sum\limits_{j}{a_{j}\sigma_{j}\frac{1}{{p - x_{j}}}}} + {\frac{a_{j}\mu_{j}}{2ɛ}\left( {\frac{1}{{p - \left( {x_{j} - {ɛ\; {\overset{\rightarrow}{n}}_{j}}} \right)}} - \frac{1}{{p - \left( {x_{j} + {ɛ\; {\overset{\rightarrow}{n}}_{j}}} \right)}}} \right)}},} \right.}} & (31)\end{matrix}$

where the source panel is converted to a point with area a_(i) and thecoefficients y_(lm) and b_(ln) are a special case of Eq. (24) withr_(j)=0. For the source located at x_(j):

μ_(lmn)=a_(j)σ_(j)y_(lm)b_(ln)  (32)

The doublet produces two sources, located at x_(j)−ϵ{right arrow over(n)}_(j) and at x_(j)+ϵ{right arrow over (n)}_(j), with farcoefficients:

$\begin{matrix}{{\mu_{lmn} = {\frac{a_{j}\mu_{j}}{2ɛ}y_{lm}b_{\ln}}}{\mu_{lmn} = {{- \frac{a_{j}\mu_{j}}{2ɛ}}y_{lm}{b_{\ln}.}}}} & (33)\end{matrix}$

The coefficients y_(lm) and b_(ln) are also the ones of Eq. (24). If thesurface is undersampled, the discretization error can manifest itselfvisually as weak currents leaking in or out of boundaries, usually farfrom the control points and near panel edges. This issue can be reducedfor inward leakage by imposing a no-through condition per point, and foroutward leakage by using the maximum principle: the magnitude of theharmonic field in the entire space is bound by the harmonic field'svalue on the surface. Near the surface, it may be assumed that

{right arrow over (n)} _(i) ·{right arrow over (u)}(p)<=argmax({rightarrow over (n)} _(i)·({right arrow over (s)}(p _(i))−{right arrow over(v)}(p _(i)))).

D. Stretching

While the previous sections focus on computing ((∇×)⁻¹{right arrow over(ω)}, the following sections focus on the dynamic terms in theright-hand side of Eq. (2). The term ({right arrow over (ω)}·∇){rightarrow over (v)} in Eq. (2) models the stretching of vorticity, and isresponsible for changing the vorticle orientation as well as introducingincreasingly high frequency velocities by transferring large scaleeddies to smaller scale eddies [Frisch and Kolmogorov 1995]. This termdoes not have an explicit analogue in Eq. (1), and comes from theexpansion into partial derivatives of the acceleration's curl. It is thechange of vorticity by the velocity gradient. Stretching may be measuredassuming that {right arrow over (w)}_(i) and {right arrow over(ω)}(p_(i)) have the same direction, which is true for an isolatedvorticle, or when r_(i)→0. As the approach taken by Zhang and Bridson[2014], two points may be displaced on either side of p_(i). Reducingthis measurement to two points may be a valid approach because the flowis incompressible, and it can be safely assumed that a stretch along{right arrow over (w)}_(i) is accompanied by a corresponding squash onthe plane perpendicular to {right arrow over (w)}_(i). After one timestep, {right arrow over (ω)}_(i)(p_(i)) becomes {right arrow over(ω)}_(i′)(p_(i′)):

$\begin{matrix}{{{{\overset{\rightarrow}{\omega}}_{i}\left( p_{i} \right)} = {2r_{i}^{{- 3}/2}{\overset{\rightarrow}{w}}_{i}}}{{{\overset{\rightarrow}{\omega}}_{i}^{\prime}\left( p_{i}^{\prime} \right)} = {2{r_{i}^{{- 3}/2}\left( {{\overset{\rightarrow}{w}}_{i} + {\frac{\Delta_{t}{{\overset{\rightarrow}{w}}_{i}}}{r_{i}}\left( {{\overset{\rightarrow}{u}\left( q_{1} \right)} - {\overset{\rightarrow}{u}\left( q_{0} \right)}} \right)}} \right)}}}{q_{0} = {p_{i} - \frac{r_{i}{\overset{\rightarrow}{w}}_{i}}{2{{\overset{\rightarrow}{w}}_{i}}}}}{q_{1} = {p_{i} + \frac{r_{i}{\overset{\rightarrow}{w}}_{i}}{2{{\overset{\rightarrow}{w}}_{i}}}}}} & (34)\end{matrix}$

To apply this measure of stretching to a vorticle, closed-form formulasmay be developed for the ellipsoidal fields of a vorticle under uniformstretching:

{right arrow over (v)} _(i)({right arrow over (x)}, s)=sξ _(i)(∥σ({rightarrow over (x)}, s)∥³ {right arrow over (w)} _(i) ×{right arrow over(x)}

{right arrow over (ω)}_(i)({right arrow over (x)}, s)=sξ _(i)(∥σ({rightarrow over (x)}, s)∥)³(2{right arrow over (w)} _(i)−3ξ_(i)(∥σ({rightarrow over (x)}, s)∥)²σ({right arrow over (x)}, s ²)×({right arrow over(w)} _(i) ×{right arrow over (x)})),  (35)

where

${\sigma \left( {\overset{\rightarrow}{x},s} \right)} = {\frac{1}{{{\overset{\rightarrow}{w}}_{i}}^{2}}{\left( {{{s^{{- 4}/5}\left( {\overset{\rightarrow}{x} \cdot {\overset{\rightarrow}{w}}_{i}} \right)}{\overset{\rightarrow}{w}}_{i}} + {{s^{7/10}\left( {{\overset{\rightarrow}{w}}_{i} \times \overset{\rightarrow}{x}} \right)} \times {\overset{\rightarrow}{w}}_{i}}} \right).}}$

When the stretching factor increases, it stretches {right arrow over(ω)}_(i) by a factor s along {right arrow over (w)}_(i), and squashes{right arrow over (ω)}_(i) by √{square root over (1/s)} along anydirection perpendicular to {right arrow over (w)}_(i), in accordancewith the deformation induced by an incompressible flow:

$\begin{matrix}{{\frac{{{\overset{\rightarrow}{\omega}}_{i}\left( {\overset{\rightarrow}{x},s} \right)} \cdot {\overset{\rightarrow}{w}}_{i}}{{{\overset{\rightarrow}{\omega}}_{i}\left( {\sigma \left( {\overset{\rightarrow}{x},s} \right)} \right)} \cdot {\overset{\rightarrow}{w}}_{i}} = s}{\frac{{{\overset{\rightarrow}{\omega}}_{i}\left( {\overset{\rightarrow}{x},s} \right)} \cdot \left( {p \times {\overset{\rightarrow}{w}}_{i}} \right)}{{{\overset{\rightarrow}{\omega}}_{i}\left( {\sigma \left( {\overset{\rightarrow}{x},s} \right)} \right)} \cdot \left( {p \times {\overset{\rightarrow}{w}}_{i}} \right)} = {\sqrt{1/s}.}}} & (36)\end{matrix}$

The rotation part of stretching may be applied to {right arrow over(w)}_(i′) and the scale part to s_(i′), and obtain a stretched vorticle:

$\begin{matrix}{r_{i}^{\prime} = r_{i}} & (37) \\{{\overset{\rightarrow}{w}}_{i}^{\prime} = {{{\overset{\rightarrow}{w}}_{i}}\frac{{\overset{\rightarrow}{\omega}}_{i}^{\prime}\left( p_{i} \right)}{{{\overset{\rightarrow}{\omega}}_{i}^{\prime}\left( p_{i}^{\prime} \right)}}}} & (38) \\{s_{i}^{\prime} = \frac{{{\overset{\rightarrow}{\omega}}_{i}^{\prime}\left( p_{i}^{\prime} \right)}}{{{\overset{\rightarrow}{\omega}}_{i}\left( p_{i} \right)}}} & (39)\end{matrix}$

The stretchable vorticle may then be unstretched in a manner thatpreserves both E_(i) and the enstrophy Ω_(i) of the stretched vorticle.Preserving enstrophy is a local way to reduce as much as possible thedisruption of this step over the vorticity dynamics.

$\begin{matrix}{\Omega_{i} = {{\frac{1}{2}{\int_{\;}^{\;}{{\overset{\rightarrow}{\omega}}_{i}^{2}{dx}}}} = {\frac{3{\pi^{2}\left( {1 + {4s_{i}^{\prime 3}}} \right)}}{64r_{i}^{\prime 3}s_{i}^{{\prime 8}/5}}{{{\overset{\rightarrow}{w}}_{i}^{\prime}}^{2}.}}}} & (40)\end{matrix}$

The stretching factor is restored to 1, and the resulting vorticle isunstretched:

$\begin{matrix}{{r_{i}^{''} = {r_{i}{s_{i}^{{\prime 4}/5}\left( \frac{5}{1 + {4s_{i}^{\prime 3}}} \right)}^{1/2}}}{{\overset{\rightarrow}{w}}_{i}^{''} = {\sqrt{\frac{r_{i}^{''}}{r_{i}}}{{\overset{\rightarrow}{w}}_{i}}\frac{{\overset{\rightarrow}{w}}_{i}^{\prime}\left( p_{i}^{\prime} \right)}{{{\overset{\rightarrow}{w}}_{i}^{\prime}\left( p_{i}^{\prime} \right)}}}}} & (41)\end{matrix}$

This can be an important new insight: when a vorticle's stretchingfactor is greater than 1, its strength must be reduced so that the meanvorticity and enstrophy can be preserved. Limit resolutions may also bespecified with a lower threshold r_(min) and upper threshold r_(max) onthe vorticle size r_(i). This limit resolution loses enstrophy, but doesnot lose energy because of Eq. (9). Thus Eq. (34), Eq. (39), and Eq.(41) can provide a way to apply stable energy preserving stretching to avorticle. In a real fluid, turbulence approaching the fluid's Kolmogorovlength mix with an effective viscosity. When a vorticle is squeezedbelow r_(min), the stretching may be turned into a higher viscositycoefficient for that vorticle, with the diffusion per vorticle of thenext section. This is necessary to remove vorticles with smallcontribution.

FIGS. 15A-15C illustrate a first vorticle 1510 (FIG. 15A) is stretchedinto a second vorticle 1520 (FIG. 15B), where the mean energy ispreserved while the enstrophy changes; and the second vorticle 1520 isthen resampled in an unstretched third vorticle 1530 (FIG. 15C) with thesame mean energy and enstrophy as the second vorticle 1520.

E. Viscosity

The term v∇²{right arrow over (ω)} in Eq. (2) models viscosity. It isthe analogue of the viscosity term of Eq. (1), with the difference thatν is varying if the density ρ is varying. Since vorticles have a sizeand aren't singular at their center, a well behaved vorticle derivativethat is consistent with the dynamics can be computed:

∇²{right arrow over (ω)}_(i)({right arrow over (x)})=15(∥{right arrowover (x)}∥ ² ξ_(i) ²−1)ξ_(i) ⁵(2{right arrow over (w)} _(i)−7ξ_(i) ²{right arrow over (x)}×{right arrow over (w)} _(i) ×{right arrow over(x)})  (42)

This equation may be solved in a multiscale manner with a near and fardiffusion in an octree, as discussed above.

$\frac{d\; {\overset{\rightarrow}{\omega}}_{i}}{dt} = {v{\nabla^{2}{\overset{\rightarrow}{\omega}}_{i}}}$

may be solved more simply, by updating the vorticle enstrophy to matchthat of {right arrow over (ω)}_(i)+Δ_(i)ν∇^(2{right arrow over (ω)})_(i). Since vorticles have no mass, this model disregards the mass ofthe advected quantity. Its advantage however is in providing users witha viscosity that is spatially tweakable per vorticle, and bound to thedynamics. The linear approximation is only valid when the diffusionΔ_(i)ν is small enough. Therefore, the diffusion may be clamped toenforce monotonicity. After one time step, the vorticle strengthbecomes:

$\begin{matrix}{{\overset{\rightarrow}{w}}_{i} = \left\{ \begin{matrix}{{\sqrt{1 + {\frac{\Delta_{t}v}{r_{i}^{2}}\left( {\frac{11025\mspace{14mu} \Delta_{t}v}{256\mspace{14mu} r_{i}^{2}} - \frac{35}{4}} \right)}}\mspace{14mu} {\overset{\rightarrow}{w}}_{i}\mspace{14mu} {if}\mspace{14mu} \Delta_{t}v} < {\frac{32}{315}r_{i}^{2}}} \\{\frac{\sqrt{5}}{3}\mspace{14mu} {\overset{\rightarrow}{w}}_{i}\mspace{14mu} {otherwise}}\end{matrix} \right.} & (43)\end{matrix}$

Another limitation of this viscous model is that it does not carry thediffusion of direction from vorticle to vorticle.

F. Buoyancy

${{\nabla{\log (\rho)}} \times \left( {\overset{\rightarrow}{F} - \frac{d\; \overset{\rightarrow}{u}}{dt}} \right)} + {\nabla{\times \overset{\rightarrow}{F}\mspace{14mu} {in}\mspace{14mu} {{Eq}.\mspace{14mu} (2)}}}$

When p is outside of solid objects, the term

${{{\nabla{\log (\rho)}} \times \left( {\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - \frac{d\; \overset{\rightarrow}{u}}{dt}} \right)} + {\nabla{\times \overset{\rightarrow}{e}}}},$

becomes and models the vorticity induced by buoyancy. This component mayrequire a varying density in order to have any effect.

New vorticles are produced from the density field ρ releasing potentialenergy. ρ may be defined with a set of particles carrying density and anambient density ρ_(A)>0, so the total density field ρ is strictlygreater than 0, as required by Eq. (2):

ρ=ρ_(A)+Σρ_(j).  (44)

FIG. 16 illustrates a density particle emits a ring of vorticles (red).The ring may be sampled with vorticles of strengths that match the massof the density particles, as visualized on points (green lattice).

The falloff ρ_(i) of a density particle j is centered at x_(j), anddefined per Appendix 7 in local coordinates q_(j)=p−x_(j) as:

$\begin{matrix}{{\rho_{j} = {m_{j}\frac{{\exp\left( \left( {1 + {\kappa_{1}\frac{q_{j}q_{j}}{2r_{j}^{2}}}} \right)^{- \kappa_{2}} \right)} - 1}{e - 1}}},} & (45)\end{matrix}$

where m_(j) is a multiplier of the particle density field, and κ_(i) isis given in Eq. (55). The newly induced vorticles are located along aring of diameter r_(j) perpendicular to

$\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - {\frac{d\overset{\rightarrow}{u}}{dt}.}$

Instead of advecting an additional filament representation, the ringwith n new vorticles may be discretized, equidistant for simplicity, andwhere n=2 in practice. Orthogonal unit vectors {right arrow over (a)}and {right arrow over (b)} may be defined such that {right arrow over(a)}×{right arrow over (b)} has the direction of

$\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - {\frac{d\overset{\rightarrow}{u}}{dt}.}$

Given n randomly selected samples α_(j), the new vorticles are in thedensity particle's local coordinates:

$\begin{matrix}{\mspace{79mu} {{x_{\alpha \; j} = {x_{j} + {\frac{r_{j}}{2}\left( {{{\cos\left( \; \alpha_{j} \right)}\overset{\rightarrow}{a}} + {{\sin \left( \alpha_{j} \right)}\overset{\rightarrow}{b}}} \right)}}}{{\overset{\rightarrow}{w}}_{\alpha \; j} = {r_{j}\sqrt{r_{j}}\kappa_{0}\frac{\pi}{n}{{\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - \frac{d\overset{\rightarrow}{u}}{dt}}}\frac{{\rho_{j}\left( x_{\alpha} \right)} + \frac{m_{j}}{e - 1}}{\rho \left( x_{\alpha} \right)}\left( {{{- {\sin \left( \alpha_{j} \right)}}\overset{\rightarrow}{a}} + {{\cos \left( \alpha_{j} \right)}\overset{\rightarrow}{b}}} \right)}}\mspace{20mu} {r_{\alpha_{j}} = r_{j}}}} & (46)\end{matrix}$

To evaluate the acceleration

$\frac{d\overset{\rightarrow}{u}}{dt},$

the velocity may be sampled both after and before the advection of theparticle, before emission of new vorticles in order to avoid thetemporal jolt of the new vorticles.

G. Vorticity Shedding

When p is at the boundary of a solid object, the term

${\bigtriangledown \; {\log (\rho)} \times \left( {\overset{\rightarrow}{F} - \frac{d\; \overset{\rightarrow}{u}}{dt}} \right)} + {\bigtriangledown \times \overset{\rightarrow}{F}}$

is Eq. (2) measures the change of vorticity at the moving object'sboundary. This component is absent in inviscid fluids, and alwayspresent in real fluids. The viscosity coefficient for vorticity sheddingcan be specified independently from the viscosity discussed above.

The surface vorticity spreads into the flow by a diffusion proportionalto viscosity coefficient ν. We show in Appendix 8 that the surfacevorticity that satisfies our boundary conditions is

$\left( {\overset{\rightarrow}{f} - \overset{\rightarrow}{e} - \frac{d\overset{\rightarrow}{u}}{dt}} \right) \times {\overset{\rightarrow}{n}.}$

The vorticity shedding may be the solution to a differential equationwith boundary condition:

$\begin{matrix}\left\{ \begin{matrix}{\frac{d\overset{\rightarrow}{\omega}}{dt} = {v\; \bigtriangledown^{2}\overset{\rightarrow}{\omega}}} \\{\overset{\rightarrow}{\omega}\underset{p \in {\delta\Omega}}{=}{{\Delta_{t}\left( {\overset{\rightarrow}{f} - \overset{\rightarrow}{e} - \frac{d\overset{\rightarrow}{u}}{dt}} \right)} \times \overset{\rightarrow}{n}}}\end{matrix} \right. & (47)\end{matrix}$

Eq. (47) may be solved by shedding vorticles. The surface may be dividedinto n panels of area a_(i) and centroid x_(i), and emit per panel avorticle that roughly approximates the heat kernel during a time stepΔ_(i):

$\begin{matrix}{{{\overset{\rightarrow}{w}}_{i} = {a_{i}\frac{r^{1/2}}{4\; \pi}\left( {\overset{\rightarrow}{f} - \overset{\rightarrow}{e} - \frac{d\overset{\rightarrow}{u}}{dt}} \right) \times \overset{\rightarrow}{n}}}{r_{i} = {\sqrt{2}{\sqrt{\Delta_{t}v}.}}}} & (48)\end{matrix}$

${\left( {\overset{\rightarrow}{f} - \overset{\rightarrow}{e} - \frac{d\overset{\rightarrow}{u}}{dt}} \right) \times \overset{\rightarrow}{n}}$

normalized over the surface may be used, as a probability distributionfunction to create samples at the locations that most affect the flow.Note that {right arrow over (f)} is the surface acceleration, as opposedto the velocity. To compute the acceleration, Eq. (6) at the previoustime may be stored on the surface, and

$\frac{d\overset{\rightarrow}{u}}{dt} \simeq \frac{{\overset{\rightarrow}{u}(t)} - {\overset{\rightarrow}{u}\left( {t - \Delta_{t}} \right)}}{\Delta_{t}}$

may be used.

This may produce the expected behavior as illustrated in FIGS. 17A-17B.FIG. 17A illustrates a constant flow moves at 10 ms⁻¹ with ν=0.1 arounda static pole of radius 1 (red): the surface sheds vorticles (blackarrows). FIG. 17B illustrates a slice of the vorticity induced by thevorticles that reveals the emergent behavior of a von Kármán vortexstreet. Furthermore, users can modify separately shedding induced bysurface acceleration {right arrow over (f)} and fluid acceleration

$\frac{d\overset{\rightarrow}{u}}{dt}.$

Optimization

User emitters, buoyancy, and shedding introduce an increasing number ofvorticles in the simulation. To address this rise of complexity, afusing step may be used. Using the existing octree data structure,similar vorticles in each octree cell may be merged. The extreme fusioncase is one vorticle per cell. The fusing method may favor the vorticleof high energy, and equals the combined energy when {right arrow over(w)}₀={right arrow over (w)}₁ and r₀=r₁:

$p = \frac{{\sqrt{E_{0}}p_{0}} + {\sqrt{E_{1}}p_{1}}}{\sqrt{E_{0}} + \sqrt{E_{1}}}$$r = \frac{{\sqrt{E_{0}}r_{0}} + {\sqrt{E_{1}}r_{1}}}{\sqrt{E_{0}} + \sqrt{E_{1}}}$$\overset{\rightarrow}{w} = {{\left( {r_{0}/r} \right)^{3/2}{\overset{\rightarrow}{w}}_{0}} + {\left( {r_{1}/r} \right)^{3/2}{\overset{\rightarrow}{w}}_{1}}}$

H. Multi-Scale Density Advection

Modern renderers handle volumes in the form of voxel grids. It maytherefore be necessary to output the density and velocity in that form.A simple method may be used as a starting point for a full solution,compatible with available volume data structures [Museth 2013; Wrenninge2016]. According to some embodiments, given a single smoke densityemitter, the space may be split into multiple grids: density is emittedin a grid with a resolution based on the position in camera space of theemitter at the times of emission. The velocity is then evaluated in thegrid, and used for advecting the density (i.e., “moving” the density).Since all densities share a common velocity field, the dynamics mayremain consistent, and multiple density resolutions may overlap in anadditive manner.

FIGS. 18A-18E illustrate schematically a method of forming a densitygrid according to some embodiments. Referring to FIG. 18A, a camera islocated at the point p. The three-dimensional space surrounding thecamera may be split into a set of regions defined by concentric spherescentered at p. In some embodiments, the radii of the spheres mayincrease exponentially as r_(j)=2^(j)×r, where r is the radius of theinner-most sphere. A region S_(j) may be defined as the space locatedbetween two consecutive spheres, i.e., a sphere of radius 2^(j)×r andthe next larger sphere of radius 2^(j+1)×r. For instance, in the exampleillustrated in FIG. 18A, the set of regions may include a first regionS₀ defined by the inner-most sphere 1810 of radius r, a second region S₁located between the inner-most sphere 1810 and the next larger sphere1820 of radius 2r, and a third region S₂ located between the sphere 1820and the next larger sphere 1830 of radius 4r, and so on and so forth.

In some embodiments, each region S_(j) may be associated with two voxelgrids: a density voxel grid and a velocity voxel grid. The size of thedensity voxel grid for each region S_(j) may increase exponentially as2^(j)×SizeA. For instance, in the examples illustrated in FIGS. 18B-18D,the first region S₀ may be associated with a 4×4 density voxel grid witha grid size of SizeA (FIG. 18B), the second region S₁ may be associatedwith a 4×4 density voxel grid with a grid size of 2×SizeA (FIG. 18C),and the third region S₂ may be associated with a 4×4 density voxel gridwith a grid size of 4×SizeA (FIG. 18D). FIG. 18E illustrates how thedensity voxel grids of the various regions S₀, S₁, and S₂ may be tiledto form an overall sparse density voxel grid.

Similarly, the size of the velocity voxel grid for each region S_(j) mayincrease exponentially as 2^(j)×SizeB. In some embodiments, SizeB may bechosen to be greater than SizeA, so that each density voxel is coveredby a velocity voxel. For example, SizeB may be chosen to be twice ofSizeA as SizeB=2×SizeA. Thus, in the example illustrated in FIGS.18A-18E, the velocity voxel grid for each region may be a 2×2 grid. Thevelocity voxel grids of the various regions S₀, S₁, and S₂ may be tiledwith each other to form an overall sparse velocity voxel grid.

In some embodiments, as a density emitter overlaps a region, anyinactive density voxels in the region that overlap with the emitter maybe activated, and densities are added to those density voxels. Forinstance, in the example illustrated in FIG. 18E, the density emitter1840 overlaps with a first density voxel 1850 and a second density voxel1860 in the second region S₁, and a third density voxel 1870 in thethird region S₂. If the first density voxel 1850, the second densityvoxel 1860, and the third density voxel 1870 are inactive (i.e., havenot been activated in a previous frame), they will be activated anddensities are added to each one of those density voxels according to thedensity emitter 1840.

Because the resolution of the density voxel grids may vary from oneregion to another, two neighboring density voxel grids can haveoverlapping density voxels. In some embodiments, to find the density ata point in a region with overlapping density voxels, the densities ofthe overlapping density voxels may be summed by the renderer.

In some embodiments, as a density voxel is activated, a correspondingvelocity voxel that covers the density voxel may be activated.Therefore, there is always active velocity voxels covering activedensity voxels. Each activated velocity voxel is then initialized bycomputing the velocity induced by the vorticles, using the accelerationoctree structure, as discussed above. The velocity voxel grid is thenused to advect the density voxel grid to the next frame, using asemi-Lagrangian advection. The term “advect” may refer to “moving”(e.g., the density).

In some embodiments, one or more density voxels may be pre-activated fora next frame using the velocity voxel grid. For example, a density voxelmay be deformed using the velocity voxel grid to predict a deformedshape of the density emitter for the next frame. Those inactive densityvoxels that overlap with the deformed density emitter may bepre-activated.

In some other embodiments, the velocity voxel grid may be optional. Somerenders may be capable of rendering sub-frame density, that is,sub-frame density voxel grids that can capture the variation over thecourse of a single frame.

I. Methods of Computer Graphic Simulation of Smoke

According to some embodiments, methods of simulating smoke in a computergraphic system may include the following steps:

(a) Add emitted vorticles, e.g., by splitting the emitter's energy amongvorticles using Eq. (10).(b) Build the octree of vorticles: store each vorticle in the cell ofsize

$2^{\lceil\frac{{logr}_{i}}{\log \; 2}\rceil},$

and allocate the far interacting cells, as described above.(c) Compute the advected field near coefficients, as described above.(d) Evaluate the advected field {right arrow over (v)} at every panellocation, using Eq. (12), as described above.(e) Compute the source and doublet panel coefficients σ_(i) and μ_(i) bysolving Eq. (27), as described above.(f) Build the octree of panels: store each panel in the cell of size

$2^{\lceil\frac{{logr}_{i}}{\log \; 2}\rceil},$

and allocate the far interacting cells, as described above.(g) Compute the harmonic field near coefficients, as described above.(h) Insert new density in the density voxels.(i) Allocate velocity voxels and compute velocity field a at everyvoxel.(j) Advect the density using the velocity voxels.(k) Compute the velocity field a at every vorticle. Note that this stepdoes not use the velocity voxels, but instead the octree structure.(l) Advect the vorticles, as described above.(m) Stretch vorticles, as described above.(n) Diffuse vorticles, as described above.

FIG. 19 is a flowchart illustrating a method of computer graphic imageprocessing for generating rendered image frames of incompressible gascorresponding to a camera view of a three-dimensional virtual spaceaccording to some embodiments.

At 1902, a plurality of vorticles is generated. Each respective vorticleis centered at a respective position in the virtual space and has arespective spatial size.

At 1904, advected field of the plurality of vorticles is computed bysolving a vorticity equation using a fast multipole method (FMM).

At 1906, the virtual space is divided into a plurality of concentricregions surrounding a position of a camera. In some embodiments, theplurality of concentric regions may be defined by a plurality ofconcentric spheres, and each respective region is located between twoconsecutive spheres, as discussed above in relation to FIG. 18A. In someother embodiments, the plurality of concentric regions may be defined bya plurality of concentric cubes, rectangular shapes, or the like.

At 1908, a respective density voxel grid is defined for each respectiveregion of the plurality of concentric regions. Each respective densityvoxel grid may have a respective first voxel size that increases withincreasing distance from the position of the camera. In someembodiments, the first voxel size may increase with increasing distancefrom the position of the camera as an exponential function, as discussedabove in relation to FIGS. 18B-18E. In other embodiments, the firstvoxel size may be constant in various regions.

At 1910, a respective velocity voxel grid is defined for each respectiveregion of the plurality of concentric regions. Each respective velocityvoxel grid may have a respective second voxel size that increases withincreasing distance from the position of the camera, as discussed above.In some embodiments, the respective second voxel size of each respectivevelocity voxel grid is greater than the respective first voxel size ofthe corresponding density voxel grid, as discussed above.

At 1912, a new density representing a smoke emitter is added. the newdensity has a first distribution in a first region of the virtual space.The first region may overlap with one or more regions of the pluralityof concentric regions.

At 1914, each density voxel of the one or more regions that overlapswith the first region is activated, and adding density is added thereto.

At 1916, each velocity voxel of the one or more regions that overlapswith the active density voxels is activated, and velocity field at thevelocity voxel is computed.

At 1918, the new density is moved based on the velocity field of eachactive velocity voxel to obtain a second distribution of the new densityin a second region of the virtual space for a next frame.

At 1920, each density voxel that overlaps with the second region isactivated, and density is added thereto.

At 1922, the density of each active density voxel and the velocity fieldof each active velocity voxel is stored in a computer-readable medium asthe next frame for display at a screen.

In some embodiments, computing the advected field of the plurality ofvorticles using the FMM may include building an octree of vorticles. Theoctree may include a plurality octree levels. Each respective octreelevel may have a cell size that is proportional to 2^(−L), where L is alevel number of the respective octree level. Each respective vorticle ofthe plurality of vorticles may be stored in a respective cell of arespective octree level having a cell size that is proportional to thespatial size of the respective vorticle. For example, each vorticle maybe stored in a cell having a cell size of of

$2^{\lceil\frac{{logr}_{i}}{\log \; 2}\rceil},$

where r_(i) is the spatial size of the respective vorticle.

A first cell containing the position of the camera and cells that aredirect neighbors of the first cell may be referred to as near cells, andall other cells may be referred to as far cells. The advected field maybe computed using the FMM by summing the vorticles stored in eachrespective far cell.

In some embodiments, velocity field of the plurality of vorticles may becomputed, and the plurality of vorticles may be moved to a next framebased on the velocity field using the octree.

J. Examples

FIGS. 20A-20B illustrate a smoke column simulated according to someembodiments. FIG. 20A shows the vorticles, and FIG. 20B shows densityadvected using the vorticles. As illustrated, turbulence size can beexplicit. For example, the turbulence may become naturally smaller asthe fluid evolves. This can also be controlled artificially. Also, thesize of the turbulence can vary from large to small. In fact, thesimulation methods disclosed herein can be multiscale and can combinetiny vorticles with very large vorticles in a stable manner. This can beadvantageous for sampling fluids in a camera frustum.

FIGS. 21A-21D illustrate an example of a simulator that may takeadvantage of a camera frustum according to some embodiments. As thesmoke emitter moves away from the camera (as shown in FIG. 21A), thereare fewer vorticles (as shown in FIG. 21B), and the number of densityvoxels in the distance may be reduced (as shown in FIG. 21C). FIG. 21Dshows a top down view of FIG. 21C, where one may more clearly see thatvoxels near the camera are smaller than voxels far from the camera.

FIG. 22A illustrates how tiny vorticles (shown in white) may interactwith large vorticles (shown in black) according to some embodiments.FIG. 22B shows the density advected using the vorticles illustrated inFIG. 22A.

FIGS. 23A-23B illustrate how the harmonic field may play the same roleas the pressure term in the classic grid-based Navier-Stokes equationaccording to some embodiments. The harmonic field may ensure that theflow does not cross boundaries of an object. The harmonic field may becomputed using the source and doublet panel method as discussed above.

As described above, embodiments of the present invention provide methodsof solving the Navier-Stokes' differential equation using sphericalharmonics across multiple resolutions. There is no constraint on placingvorticles of arbitrary size. The methods can be very stable with mixedscales. The divergence-free basis functions may produce a field with nodivergence by construction. The sampling resolution of density, dynamicsand boundary condition may be decoupled from each other, and each can besampled adaptively. The panel method may be further improved byfiltering the advected field with a filter size proportional to thepanel size. Vorticles that carry a volume of vorticity can be advected,as opposed to regularizing a point vorticity. This may lead to stableand energy preserving stretching. The stretching model may be furtherimproved by splitting stretched vorticles along the stretchingdirections, which would be closer to the physical phenomenon thatspreads turbulence along stretched features. It may also be advantageousto vary the multipole expansion across scale for better control of theerror. Additional stability may be achieved by integrating twistsinstead of velocity. Also, the model for buoyancy, which may play a keyrole in gas motion, is proposed for producing more lively smokes andpyroclastic effects. The methods may provide the full set of featuresexpected from a gas simulation, and can be an important tool for makingvisual effects in computer graphic simulation.

K. Computer Systems

FIG. 24 is a simplified block diagram of system 2400 for creatingcomputer graphics imagery (CGI) and computer-aided animation that mayimplement or incorporate various embodiments. In this example, system2400 can include one or more design computers 2410, object library 2420,one or more object modeler systems 2430, one or more object articulationsystems 2440, one or more object animation systems 2450, one or moreobject simulation systems 2460, and one or more object rendering systems2470. Any of the systems 2430-870 may be invoked by or used directly bya user of the one or more design computers 2410 and/or automaticallyinvoked by or used by one or more processes associated with the one ormore design computers 2410. Any of the elements of system 2400 caninclude hardware and/or software elements configured for specificfunctions.

The one or more design computers 2410 can include hardware and softwareelements configured for designing CGI and assisting with computer-aidedanimation. Each of the one or more design computers 2410 may be embodiedas a single computing device or a set of one or more computing devices.Some examples of computing devices are PCs, laptops, workstations,mainframes, cluster computing system, grid computing systems, cloudcomputing systems, embedded devices, computer graphics devices, gamingdevices and consoles, consumer electronic devices having programmableprocessors, or the like. The one or more design computers 2410 may beused at various stages of a production process (e.g., pre-production,designing, creating, editing, simulating, animating, rendering,post-production, etc.) to produce images, image sequences, motionpictures, video, audio, or associated effects related to CGI andanimation.

In one example, a user of the one or more design computers 2410 actingas a modeler may employ one or more systems or tools to design, create,or modify objects within a computer-generated scene. The modeler may usemodeling software to sculpt and refine a neutral 3D model to fitpredefined aesthetic needs of one or more character designers. Themodeler may design and maintain a modeling topology conducive to astoryboarded range of deformations. In another example, a user of theone or more design computers 2410 acting as an articulator may employone or more systems or tools to design, create, or modify controls oranimation variables (avars) of models. In general, rigging is a processof giving an object, such as a character model, controls for movement,therein “articulating” its ranges of motion. The articulator may workclosely with one or more animators in rig building to provide and refinean articulation of the full range of expressions and body movementneeded to support a character's acting range in an animation. In afurther example, a user of design computer 2410 acting as an animatormay employ one or more systems or tools to specify motion and positionof one or more objects over time to produce an animation.

Object library 2420 can include elements configured for storing andaccessing information related to objects used by the one or more designcomputers 2410 during the various stages of a production process toproduce CGI and animation. Some examples of object library 2420 caninclude a file, a database, or other storage devices and mechanisms.Object library 2420 may be locally accessible to the one or more designcomputers 2410 or hosted by one or more external computer systems.

Some examples of information stored in object library 2420 can includean object itself, metadata, object geometry, object topology, rigging,control data, animation data, animation cues, simulation data, texturedata, lighting data, shader code, or the like. An object stored inobject library 2420 can include any entity that has an n-dimensional(e.g., 2D or 3D) surface geometry. The shape of the object can include aset of points or locations in space (e.g., object space) that make upthe object's surface. Topology of an object can include the connectivityof the surface of the object (e.g., the genus or number of holes in anobject) or the vertex/edge/face connectivity of an object.

The one or more object modeling systems 2430 can include hardware and/orsoftware elements configured for modeling one or more objects. Modelingcan include the creating, sculpting, and editing of an object. Invarious embodiments, the one or more object modeling systems 2430 may beconfigured to generated a model to include a description of the shape ofan object. The one or more object modeling systems 2430 can beconfigured to facilitate the creation and/or editing of features, suchas non-uniform rational B-splines or NURBS, polygons and subdivisionsurfaces (or SubDivs), that may be used to describe the shape of anobject. In general, polygons are a widely used model medium due to theirrelative stability and functionality. Polygons can also act as thebridge between NURBS and SubDivs. NURBS are used mainly for theirready-smooth appearance and generally respond well to deformations.SubDivs are a combination of both NURBS and polygons representing asmooth surface via the specification of a coarser piecewise linearpolygon mesh. A single object may have several different models thatdescribe its shape.

The one or more object modeling systems 2430 may further generate modeldata (e.g., 2D and 3D model data) for use by other elements of system2400 or that can be stored in object library 2420. The one or moreobject modeling systems 2430 may be configured to allow a user toassociate additional information, metadata, color, lighting, rigging,controls, or the like, with all or a portion of the generated modeldata.

The one or more object articulation systems 2440 can include hardwareand/or software elements configured to articulating one or morecomputer-generated objects. Articulation can include the building orcreation of rigs, the rigging of an object, and the editing of rigging.In various embodiments, the one or more articulation systems 2440 can beconfigured to enable the specification of rigging for an object, such asfor internal skeletal structures or eternal features, and to define howinput motion deforms the object. One technique is called “skeletalanimation,” in which a character can be represented in at least twoparts: a surface representation used to draw the character (called theskin) and a hierarchical set of bones used for animation (called theskeleton).

The one or more object articulation systems 2440 may further generatearticulation data (e.g., data associated with controls or animationsvariables) for use by other elements of system 2400 or that can bestored in object library 2420. The one or more object articulationsystems 2440 may be configured to allow a user to associate additionalinformation, metadata, color, lighting, rigging, controls, or the like,with all or a portion of the generated articulation data.

The one or more object animation systems 2450 can include hardwareand/or software elements configured for animating one or morecomputer-generated objects. Animation can include the specification ofmotion and position of an object over time. The one or more objectanimation systems 2450 may be invoked by or used directly by a user ofthe one or more design computers 2410 and/or automatically invoked by orused by one or more processes associated with the one or more designcomputers 2410.

In various embodiments, the one or more animation systems 2450 may beconfigured to enable users to manipulate controls or animation variablesor utilized character rigging to specify one or more key frames ofanimation sequence. The one or more animation systems 2450 generateintermediary frames based on the one or more key frames. In someembodiments, the one or more animation systems 2450 may be configured toenable users to specify animation cues, paths, or the like according toone or more predefined sequences. The one or more animation systems 2450generate frames of the animation based on the animation cues or paths.In further embodiments, the one or more animation systems 2450 may beconfigured to enable users to define animations using one or moreanimation languages, morphs, deformations, or the like.

The one or more object animations systems 2450 may further generateanimation data (e.g., inputs associated with controls or animationsvariables) for use by other elements of system 2400 or that can bestored in object library 2420. The one or more object animations systems2450 may be configured to allow a user to associate additionalinformation, metadata, color, lighting, rigging, controls, or the like,with all or a portion of the generated animation data.

The one or more object simulation systems 2460 can include hardwareand/or software elements configured for simulating one or morecomputer-generated objects. Simulation can include determining motionand position of an object over time in response to one or more simulatedforces or conditions. The one or more object simulation systems 2460 maybe invoked by or used directly by a user of the one or more designcomputers 2410 and/or automatically invoked by or used by one or moreprocesses associated with the one or more design computers 2410.

In various embodiments, the one or more object simulation systems 2460may be configured to enables users to create, define, or edit simulationengines, such as a physics engine or physics processing unit (PPU/GPGPU)using one or more physically-based numerical techniques. In general, aphysics engine can include a computer program that simulates one or morephysics models (e.g., a Newtonian physics model), using variables suchas mass, velocity, friction, wind resistance, or the like. The physicsengine may simulate and predict effects under different conditions thatwould approximate what happens to an object according to the physicsmodel. The one or more object simulation systems 2460 may be used tosimulate the behavior of objects, such as hair, fur, and cloth, inresponse to a physics model and/or animation of one or more charactersand objects within a computer-generated scene.

The one or more object simulation systems 2460 may further generatesimulation data (e.g., motion and position of an object over time) foruse by other elements of system 2400 or that can be stored in objectlibrary 2420. The generated simulation data may be combined with or usedin addition to animation data generated by the one or more objectanimation systems 2450. The one or more object simulation systems 2460may be configured to allow a user to associate additional information,metadata, color, lighting, rigging, controls, or the like, with all or aportion of the generated simulation data.

The one or more object rendering systems 2470 can include hardwareand/or software element configured for “rendering” or generating one ormore images of one or more computer-generated objects. “Rendering” caninclude generating an image from a model based on information such asgeometry, viewpoint, texture, lighting, and shading information. The oneor more object rendering systems 2470 may be invoked by or used directlyby a user of the one or more design computers 2410 and/or automaticallyinvoked by or used by one or more processes associated with the one ormore design computers 2410. One example of a software program embodiedas the one or more object rendering systems 2470 can includePhotoRealistic RenderMan, or PRMan, produced by Pixar Animations Studiosof Emeryville, Calif.

In various embodiments, the one or more object rendering systems 2470can be configured to render one or more objects to produce one or morecomputer-generated images or a set of images over time that provide ananimation. The one or more object rendering systems 2470 may generatedigital images or raster graphics images.

In various embodiments, a rendered image can be understood in terms of anumber of visible features. Some examples of visible features that maybe considered by the one or more object rendering systems 2470 mayinclude shading (e.g., techniques relating to how the color andbrightness of a surface varies with lighting), texture-mapping (e.g.,techniques relating to applying detail information to surfaces orobjects using maps), bump-mapping (e.g., techniques relating tosimulating small-scale bumpiness on surfaces), fogging/participatingmedium (e.g., techniques relating to how light dims when passing throughnon-clear atmosphere or air) shadows (e.g., techniques relating toeffects of obstructing light), soft shadows (e.g., techniques relatingto varying darkness caused by partially obscured light sources),reflection (e.g., techniques relating to mirror-like or highly glossyreflection), transparency or opacity (e.g., techniques relating to sharptransmissions of light through solid objects), translucency (e.g.,techniques relating to highly scattered transmissions of light throughsolid objects), refraction (e.g., techniques relating to bending oflight associated with transparency), diffraction (e.g., techniquesrelating to bending, spreading and interference of light passing by anobject or aperture that disrupts the ray), indirect illumination (e.g.,techniques relating to surfaces illuminated by light reflected off othersurfaces, rather than directly from a light source, also known as globalillumination), caustics (e.g., a form of indirect illumination withtechniques relating to reflections of light off a shiny object, orfocusing of light through a transparent object, to produce brighthighlights on another object), depth of field (e.g., techniques relatingto how objects appear blurry or out of focus when too far in front of orbehind the object in focus), motion blur (e.g., techniques relating tohow objects appear blurry due to high-speed motion, or the motion of thecamera), non-photorealistic rendering (e.g., techniques relating torendering of scenes in an artistic style, intended to look like apainting or drawing), or the like.

The one or more object rendering systems 2470 may further render images(e.g., motion and position of an object over time) for use by otherelements of system 2400 or that can be stored in object library 2420.The one or more object rendering systems 2470 may be configured to allowa user to associate additional information or metadata with all or aportion of the rendered image.

FIG. 25 is a block diagram of computer system 2500. FIG. 25 is merelyillustrative. In some embodiments, a computer system includes a singlecomputer apparatus, where the subsystems can be the components of thecomputer apparatus. In other embodiments, a computer system can includemultiple computer apparatuses, each being a subsystem, with internalcomponents. Computer system 2500 and any of its components or subsystemscan include hardware and/or software elements configured for performingmethods described herein.

Computer system 2500 may include familiar computer components, such asone or more one or more data processors or central processing units(CPUs) 2505, one or more graphics processors or graphical processingunits (GPUs) 2510, memory subsystem 2515, storage subsystem 2520, one ormore input/output (I/O) interfaces 2525, communications interface 2530,or the like. Computer system 2500 can include system bus 2535interconnecting the above components and providing functionality, suchconnectivity and inter-device communication

The one or more data processors or central processing units (CPUs) 2505can execute logic or program code or for providing application-specificfunctionality. Some examples of CPU(s) 2505 can include one or moremicroprocessors (e.g., single core and multi-core) or micro-controllers,one or more field-gate programmable arrays (FPGAs), andapplication-specific integrated circuits (ASICs). As user herein, aprocessor includes a multi-core processor on a same integrated chip, ormultiple processing units on a single circuit board or networked.

The one or more graphics processor or graphical processing units (GPUs)2510 can execute logic or program code associated with graphics or forproviding graphics-specific functionality. GPUs 2510 may include anyconventional graphics processing unit, such as those provided byconventional video cards. In various embodiments, GPUs 2510 may includeone or more vector or parallel processing units. These GPUs may be userprogrammable, and include hardware elements for encoding/decodingspecific types of data (e.g., video data) or for accelerating 2D or 3Ddrawing operations, texturing operations, shading operations, or thelike.

The one or more graphics processors or graphical processing units (GPUs)2510 may include any number of registers, logic units, arithmetic units,caches, memory interfaces, or the like.

Memory subsystem 2515 can store information, e.g., usingmachine-readable articles, information storage devices, orcomputer-readable storage media. Some examples can include random accessmemories (RAM), read-only-memories (ROMS), volatile memories,non-volatile memories, and other semiconductor memories. Memorysubsystem 2515 can include data and program code 2540.

Storage subsystem 2520 can also store information using machine-readablearticles, information storage devices, or computer-readable storagemedia. Storage subsystem 2520 may store information using storage media2545. Some examples of storage media 2545 used by storage subsystem 2520can include floppy disks, hard disks, optical storage media such asCD-ROMS, DVDs and bar codes, removable storage devices, networkedstorage devices, or the like. In some embodiments, all or part of dataand program code 2540 may be stored using storage subsystem 2520.

The one or more input/output (I/O) interfaces 2525 can perform I/Ooperations. One or more input devices 2550 and/or one or more outputdevices 2555 may be communicatively coupled to the one or more I/Ointerfaces 2525. The one or more input devices 2550 can receiveinformation from one or more sources for computer system 2500. Someexamples of the one or more input devices 2550 may include a computermouse, a trackball, a track pad, a joystick, a wireless remote, adrawing tablet, a voice command system, an eye tracking system, externalstorage systems, a monitor appropriately configured as a touch screen, acommunications interface appropriately configured as a transceiver, orthe like. In various embodiments, the one or more input devices 2550 mayallow a user of computer system 2500 to interact with one or morenon-graphical or graphical user interfaces to enter a comment, selectobjects, icons, text, user interface widgets, or other user interfaceelements that appear on a monitor/display device via a command, a clickof a button, or the like.

The one or more output devices 2555 can output information to one ormore destinations for computer system 2500. Some examples of the one ormore output devices 2555 can include a printer, a fax, a feedback devicefor a mouse or joystick, external storage systems, a monitor or otherdisplay device, a communications interface appropriately configured as atransceiver, or the like. The one or more output devices 2555 may allowa user of computer system 2500 to view objects, icons, text, userinterface widgets, or other user interface elements. A display device ormonitor may be used with computer system 2500 and can include hardwareand/or software elements configured for displaying information.

Communications interface 2530 can perform communications operations,including sending and receiving data. Some examples of communicationsinterface 2530 may include a network communications interface (e.g.Ethernet, Wi-Fi, etc.). For example, communications interface 2530 maybe coupled to communications network/external bus 2560, such as acomputer network, a USB hub, or the like. A computer system can includea plurality of the same components or subsystems, e.g., connectedtogether by communications interface 2530 or by an internal interface.In some embodiments, computer systems, subsystem, or apparatuses cancommunicate over a network. In such instances, one computer can beconsidered a client and another computer a server, where each can bepart of a same computer system. A client and a server can each includemultiple systems, subsystems, or components.

Computer system 2500 may also include one or more applications (e.g.,software components or functions) to be executed by a processor toexecute, perform, or otherwise implement techniques disclosed herein.These applications may be embodied as data and program code 2540.Additionally, computer programs, executable computer code,human-readable source code, shader code, rendering engines, or the like,and data, such as image files, models including geometrical descriptionsof objects, ordered geometric descriptions of objects, proceduraldescriptions of models, scene descriptor files, or the like, may bestored in memory subsystem 2515 and/or storage subsystem 2520.

Such programs may also be encoded and transmitted using carrier signalsadapted for transmission via wired, optical, and/or wireless networksconforming to a variety of protocols, including the Internet. As such, acomputer readable medium according to an embodiment of the presentinvention may be created using a data signal encoded with such programs.Computer readable media encoded with the program code may be packagedwith a compatible device or provided separately from other devices(e.g., via Internet download). Any such computer readable medium mayreside on or within a single computer product (e.g. a hard drive, a CD,or an entire computer system), and may be present on or within differentcomputer products within a system or network. A computer system mayinclude a monitor, printer, or other suitable display for providing anyof the results mentioned herein to a user.

Any of the methods described herein may be totally or partiallyperformed with a computer system including one or more processors, whichcan be configured to perform the steps. Thus, embodiments can bedirected to computer systems configured to perform the steps of any ofthe methods described herein, potentially with different componentsperforming a respective steps or a respective group of steps. Althoughpresented as numbered steps, steps of methods herein can be performed ata same time or in a different order. Additionally, portions of thesesteps may be used with portions of other steps from other methods. Also,all or portions of a step may be optional. Additionally, any of thesteps of any of the methods can be performed with modules, circuits, orother means for performing these steps.

The specific details of particular embodiments may be combined in anysuitable manner without departing from the spirit and scope ofembodiments of the invention. However, other embodiments of theinvention may be directed to specific embodiments relating to eachindividual aspect, or specific combinations of these individual aspects.

The above description of exemplary embodiments of the invention has beenpresented for the purposes of illustration and description. It is notintended to be exhaustive or to limit the invention to the precise formdescribed, and many modifications and variations are possible in lightof the teaching above. The embodiments were chosen and described inorder to best explain the principles of the invention and its practicalapplications to thereby enable others skilled in the art to best utilizethe invention in various embodiments and with various modifications asare suited to the particular use contemplated.

A recitation of “a”, “an” or “the” is intended to mean “one or more”unless specifically indicated to the contrary.

All patents, patent applications, publications, and descriptionsmentioned here are incorporated by reference in their entirety for allpurposes. None is admitted to be prior art.

L. References

Alexis Angelidis, Marie-Paule Cani, Geoff Wyvill, and Scott King. 2006a.Swirlingsweepers: constant volume modeling. Graphical Models (GMOD) 68,4 (July 2006). Special issue on PG'04.

Alexis Angelidis, Fabrice Neyret, Karan Singh, and Derek Nowrouzezahrai.2006b. A Controllable, Fast and Stable Basis for Vortex Based SmokeSimulation. In Proceedings of the 2006 ACM SIGGRAPH/EurographicsSymposium on Computer Animation. Eurographics Association, 25-32.

Rutherford Aris. 2012. Vectors, Tensors and the Basic Equations of FluidMechanics. Dover Publications.

J. Thomas Beale and Andrew Majda. 1982. Vortex methods I: Convergence inthree dimensions. Math. Comp. 39 (July 1982), 1-27.

Tyson Brochu, Todd Keeler, and Robert Bridson. 2012. Linear-time SmokeAnimation with Vortex Sheet Meshes. In Proceedings of the ACMSIGGRAPH/Eurographics Symposium on Computer Animation. EurographicsAssociation, 87-95.

Michael J. Carley. 2012. Analytical Formulae for Potential Integrals onTriangles. Journal of Applied Mechanics 80, 4 (October 2012).

Albert Chern, Felix Knoppel, Ulrich Pinkall, Peter Schroder, and SteffenWeißmann. 2016. Schrödinger's Smoke. ACM Trans. Graph. 35, 4, Article 77(July 2016), 13 pages.

Georges-Henri. Cottet and Petros D. Koumoutsakos. 2000. Vortex Methods:Theory and Practice. Cambridge University Press.

Benoit Couet, Oscar Buneman, and Anthony Leonard. 1981. Simulation ofthree-dimensional incompressible flows with a vortex-in-cell method. J.Comput. Phys. 39 (February 1981), 305-328.

Tyler de Witt, Christian Lessig, and Eugene Fiume. 2012. FluidSimulation Using Laplacian Eigenfunctions. ACM Trans. Graph. 31, 1,Article 10 (February 2012), 11 pages.

Sharif Elcott, Yiying Tong, Eva Kanso, Peter Schröder, and MathieuDesbrun. 2007. Stable, Circulation-preserving, Simplicial Fluids. ACMTrans. Graph. 26, 1, Article 4 (January 2007), 12 pages.

Larry L. Erickson. 1990. Panel methods: An introduction. TechnicalReport. NASA Ames Research Center.

Nick Foster and Dimitris Metaxas. 1997. Modeling the Motion of a Hot,Turbulent Gas. In Proceedings of the 24th Annual Conference on ComputerGraphics and Interactive Techniques (SIGGRAPH '97). ACMPress/Addison-Wesley Publishing Co., 181-188.

Uriel Frisch and A. N. Kolmogorov. 1995. Turbulence: The Legacy ofAndrei N. Kolmogorov. Cambridge University Press.

Robert A. Gingold and Joe J. Monaghan. 1977. Smoothed particlehydrodynamics—Theory and application to non-spherical stars. MonthlyNotices of the Royal Astronomical Society 181 (November 1977), 375-389.

Leslie Greengard and Vladimir Rokhlin. 1987. A Fast Algorithm forParticle Simulations. J. Comput. Phys. 73, 2 (December 1987), 325-348.

Doyub Kim, Seung Woo Lee, Oh-Young Song, and Ko Hyeong-Seok. 2012.Baroclinic Turbulence with Varying Density and Temperature. IEEETransactions on Visualization and Computer Graphics 18 (2012),1488-1495.

Theodore Kim, Nils Thürey, Doug James, and Markus Gross. 2008. WaveletTurbulence for Fluid Simulation. ACM Trans. Graph. 27, 3, Article 50(August 2008), 6 pages.

Ken Museth. 2013. VDB: High-resolution Sparse Volumes with DynamicTopology. ACM Trans. Graph. 32, 3 (July 2013).

Sang Il Park and Myoung Jun Kim. 2005. Vortex Fluid for GaseousPhenomena. In Proceedings of the 2005 ACM SIGGRAPH/EurographicsSymposium on Computer Animation (SCA '05). ACM, 261-270.

Tobias Pfaff, Nils Thuerey, and Markus Gross. 2012. Lagrangian VortexSheets for Animating Fluids. ACM Trans. Graph. 31, 4, Article 112 (July2012), 8 pages.

Tobias Pfaff, Nils Thuerey, Andrew Selle, and Markus Gross. 2009.Synthetic Turbulence Using Artificial Boundary Layers. ACM Trans. Graph.28, 5, Article 121 (December 2009), 10 pages.

Andrew Selle, Nick Rasmussen, and Ronald Fedkiw. 2005. A Vortex ParticleMethod for Smoke, Water and Explosions. In ACM SIGGRAPH 2005 Papers(SIGGRAPH '05). ACM, 910-914.

Tarun Kumar Sheel. 2011. Acceleration of Vortex Methods CalculationUsing FMM and MDGRAPE-3. Progress In Electromagnetics Research B 27(2011), 327-348.

SideFX. 2016. Houdini 15.5. (May 2016). https://www.sidefx.com

Jos Stam. 1999. Stable Fluids. In Proceedings of the 26th AnnualConference on Computer Graphics and Interactive Techniques (SIGGRAPH'99). ACM Press/Addison-Wesley Publishing Co., 121-128.

Adriaan van Oosterom. 2011. Closed-form analytical expressions for thepotential fields generated by triangular monolayers with linearlydistributed source strength. Medical & Biological Engineering &Computing (MBEC) 50, 1 (October 2011), 1-9.

Adriaan van Oosterom and Jan Strackee. 1983. The solid angle of a planetriangle (IEEE Transactions on Biomedial Engineering), Vol. 30. IEEE,125-126.

SteffenWeissmann and Ulrich Pinkall. 2010. Filament-based Smoke withVortex Shedding and Variational Reconnection. ACM Trans. Graph. 29, 4,Article 115 (July 2010), 12 pages.

Magnus Wrenninge. 2016. Efficient Rendering of Volumetric Motion BlurUsing Temporally Unstructured Volumes. Journal of Computer GraphicsTechniques (JCGT) 5, 1 (January 2016).

Lexing Ying. 2012. A pedestrian introduction to fast multipole methods.Science China Mathematics 55, 5 (2012), 1043-1051.

Xinxin Zhang and Robert Bridson. 2014. A PPPM Fast Summation Method forFluids and Beyond. ACM Trans. Graph. 33, 6 (November 2014).

APPENDIX 1

Nomenclature

t time Δ_(t) time step p position p_(i) vorticle center {right arrowover (u)} velocity field {right arrow over (h)} harmonic field {rightarrow over (ψ)} stream field {right arrow over (ψ)}_(i) vorticle streamfield {right arrow over (v)} advected field {right arrow over (v)}_(i)vorticle velocity field {right arrow over (ω)} vorticity field {rightarrow over (ω)}_(i) vorticle vorticity field {right arrow over (w)}_(i)vorticle strength r_(i) vorticle size s_(i) vorticle stretching factorE_(i) vorticle mean energy Ω_(i) vorticle enstrophy ρ density fieldρ_(A) ambient density field ρ_(j) particle density field m_(j) particledensity multiplier {right arrow over (s)} solid objects velocity {rightarrow over (f)} solid objects acceleration ν kinematic viscosity {rightarrow over (g)} gravity {right arrow over (e)} user external forcesY_(lm) real spherical harmonics P_(lm) Legendre polynomial {right arrowover (μ)}_(lmn) far vector coefficient {right arrow over (κ)}_(lmn) nearvector coefficient μ_(lmn) far scalar coefficient κ_(lmn) near scalarcoefficient y_(lm) spherical coefficient b_(ln) radial coefficient δΩboundary surface δΩ_(j) triangle surface σ_(j) source coefficient μ_(j)doublet coefficient S_(j) source field D_(j) doublet field

APPENDIX 2

Field Basis Functions

The first few basis functions are of the form ∥p∥^(n) (n{right arrowover (a)}_(lm)+{right arrow over (b)}_(lm)), and are given in

Cartesian coordinates:

$\mspace{20mu} {{\bigtriangledown \left( {{Y_{0,0}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{2}}\sqrt{\frac{1}{4\; \pi}}({np})}}$$\mspace{20mu} {{\bigtriangledown \left( {{Y_{1,{- 1}}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{3}}\sqrt{\frac{3}{4\; \pi}}\left( {{{- {np}_{y}}p} + {p_{y}p} - {{p}^{2}\overset{\rightarrow}{y}}} \right)}}$$\mspace{20mu} {{\bigtriangledown \left( {{Y_{1,0}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{3}}\sqrt{\frac{3}{4\; \pi}}\left( {{{np}_{z}p} - {p_{z}p} + {{p}^{2}\overset{\rightarrow}{z}}} \right)}}$$\mspace{20mu} {{\bigtriangledown \left( {{Y_{1,1}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{3}}\sqrt{\frac{3}{4\; \pi}}\left( {{{np}_{x}p} + {p_{x}p} - {{p}^{2}\overset{\rightarrow}{x}}} \right)}}$${\bigtriangledown \left( {{Y_{2,{- 2}}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{4}}\sqrt{\frac{15}{4\; \pi}}\left( {{{np}_{x}p_{y}p} + {{p}^{2}\left( {{p_{y}\overset{\rightarrow}{x}} + {p_{x}\overset{\rightarrow}{y}}} \right)} - {2\; p_{x}p_{y}p}} \right)}$${\bigtriangledown \left( {{Y_{2,{- 1}}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{4}}\sqrt{\frac{15}{4\; \pi}}\left( {{{- {np}_{y}}p_{z}p} + {2\; p_{y}p_{z}p} - {{p}^{2}\left( {{p_{z}\overset{\rightarrow}{y}} + {p_{y}\overset{\rightarrow}{z}}} \right)}} \right)}$${\bigtriangledown \left( {{Y_{2,0}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{4}}\sqrt{\frac{5}{16\; \pi}}\left( {{{- {n\left( {{3\; p_{z}^{2}} - {p}^{2}} \right)}}p} + {6\; {p_{z}\left( {{p_{z}p} - {{p}^{2}\overset{\rightarrow}{z}}} \right)}}} \right)}$${\bigtriangledown \left( {{Y_{2,1}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{4}}\sqrt{\frac{15}{4\; \pi}}\left( {{{np}_{z}p_{x}p} + {{p}^{2}\left( {{p_{x}\overset{\rightarrow}{z}} + {p_{z}\overset{\rightarrow}{x}}} \right)} - {2\; p_{z}p_{x}p}} \right)}$${\bigtriangledown \left( {{Y_{2,2}(p)}{p}^{n}} \right)} = {\frac{{p}^{n}}{{p}^{4}}\sqrt{\frac{15}{4\; \pi}}\left( {{n\frac{p_{x}^{2} - p_{y}^{2}}{2}p} + {{p}^{2}\left( {{p_{x}\overset{\rightarrow}{x}} - {p_{z}\overset{\rightarrow}{y}}} \right)} + {\left( {p_{y}^{2} - p_{x}^{2}} \right)p}} \right)}$

APPENDIX 3

First Vorticle to Far Coefficient

The far radial coefficients b_(ln) multiply the radial part of thebasis:

$\begin{matrix}{B_{l} = {\sum\limits_{n = 1}^{\infty}\; {b_{\ln}\frac{1}{r^{n}}}}} & (49)\end{matrix}$

For a cell centered a c, and d_(i)=∥p_(i)−c∥, the first few far radialcoefficients b_(ln) are:

$\begin{matrix}{{{B_{0}\left( {r,d_{i}} \right)} = {\frac{1}{r} - \frac{r_{i}^{2}}{2\; r^{3}} + \frac{{{- 4}\; d_{i}^{2}r_{i}^{2}} + {3\; r_{i}^{4}}}{8\; r^{5}} + {O\left( \frac{1}{r^{7}} \right)}}}{{B_{1}\left( {r,d_{i}} \right)} = {\frac{d_{i}}{r^{2}} - \frac{3\; d_{i}r_{i}^{2}}{2\; r^{4}} - \frac{3\; {d_{i}\left( {{4\; d_{i}^{2}r_{i}^{2}} - {5\; r_{i}^{4}}} \right)}}{8\; r^{6}} + {O\left( \frac{1}{r^{7}} \right)}}}{{B_{2}\left( {r,d_{i}} \right)} = {\frac{d_{i}^{2}}{r^{3}} - \frac{5\; d_{i}^{2}r_{i}^{2}}{2\; r^{5}} + {O\left( \frac{1}{r^{7}} \right)}}}} & (50)\end{matrix}$

The first few spherical coefficients are:

$\begin{matrix}{{y_{0,0} = \sqrt{4\; \pi}}{y_{1,{- 1}} = {{- \sqrt{\frac{4\; \pi}{3}}}\frac{p_{iy}}{p_{i}}}}{y_{1,0} = {\sqrt{\frac{4\; \pi}{3}}\frac{p_{iz}}{p_{i}}}}{y_{1,1} = {{- \sqrt{\frac{4\; \pi}{3}}}\frac{p_{ix}}{p_{i}}}}{y_{2,{- 2}} = {\sqrt{\frac{12\; \pi}{5}}\frac{p_{ix}p_{iy}}{{p_{i}}^{2}}}}{y_{2,{- 1}} = {{- \sqrt{\frac{12\; \pi}{5}}}\frac{p_{iy}p_{iz}}{{p_{i}}^{2}}}}{y_{2,0} = {{- \sqrt{\frac{\pi}{5}}}\left( {1 - {3\frac{p_{iz}^{2}}{{p_{i}}^{2}}}} \right)}}{y_{2,1} = {{- \sqrt{\frac{12\; \pi}{5}}}\frac{p_{ix}p_{iz}}{{p_{i}}^{2}}}}{y_{2,2} = {\sqrt{\frac{3\; \pi}{5}}\frac{\left( {p_{ix}^{2} - p_{iy}^{2}} \right)}{{p_{i}}^{2}}}}} & (51)\end{matrix}$

APPENDIX 4

First Far Coefficient Translation

The first few coefficients:

$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{0,0,1}^{\prime} = {\overset{\rightarrow}{\mu}}_{0,0,1}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{0,0,2}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{0,0,3}^{\prime} = {\overset{\rightarrow}{\mu}}_{0,0,3}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,{- 1},1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,{- 1},2}^{\prime} = {{\overset{\rightarrow}{\mu}}_{1,{- 1},2}^{\prime} + {\sqrt{\frac{1}{3}}c_{y}{\overset{\rightarrow}{\mu}}_{0,0,1}}}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,{- 1},3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,0,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,0,2}^{\prime} = {{\overset{\rightarrow}{\mu}}_{1,0,2}^{\prime} - {\sqrt{\frac{1}{3}}c_{z}{\overset{\rightarrow}{\mu}}_{0,0,1}}}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,0,3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,1,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,1,2}^{\prime} = {{\overset{\rightarrow}{\mu}}_{1,1,2} + {\sqrt{\frac{1}{3}}c_{x}{\overset{\rightarrow}{\mu}}_{0,0,1}}}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{1,1,3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,{- 2},1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,{- 2},2}^{\prime} = 0}$${\overset{\rightarrow}{\mu}}_{2,{- 2},3}^{\prime} = {{\overset{\rightarrow}{\mu}}_{2,{- 2},3}^{\prime} + {\sqrt{\frac{3}{5}}c_{x}c_{y}{\overset{\rightarrow}{\mu}}_{0,0,1}} + {\sqrt{\frac{9}{5}}c_{x}{\overset{\rightarrow}{\mu}}_{1,{- 1},2}} + {\sqrt{\frac{9}{5}}c_{y}{\overset{\rightarrow}{\mu}}_{1,1,2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,{- 1},1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,{- 1},2}^{\prime} = 0}$${\overset{\rightarrow}{\mu}}_{2,{- 1},3}^{\prime} = {{\overset{\rightarrow}{\mu}}_{2,{- 1},3}^{\prime} - {\sqrt{\frac{3}{5}}c_{y}c_{z}{\overset{\rightarrow}{\mu}}_{0,0,1}} - {\sqrt{\frac{9}{5}}c_{z}{\overset{\rightarrow}{\mu}}_{1,{- 1},2}} + {\sqrt{\frac{9}{5}}c_{y}{\overset{\rightarrow}{\mu}}_{1,0,2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,0,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,0,2}^{\prime} = 0}$${\overset{\rightarrow}{\mu}}_{2,0,3}^{\prime} = {{\overset{\rightarrow}{\mu}}_{2,0,3}^{\prime} - {\sqrt{\frac{1}{20}}\left( {c_{x}^{2} + c_{y}^{2} - {2\; c_{z}^{2}}} \right){\overset{\rightarrow}{\mu}}_{0,0,1}} - {\sqrt{\frac{3}{5}}c_{y}{\overset{\rightarrow}{\mu}}_{1,{- 1},2}} - {\sqrt{\frac{12}{5}}c_{z}{\overset{\rightarrow}{\mu}}_{1,0,2}} - {\sqrt{\frac{3}{5}}c_{x}{\overset{\rightarrow}{\mu}}_{1,1,2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,1,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,1,2}^{\prime} = 0}$${\overset{\rightarrow}{\mu}}_{2,1,3}^{\prime} = {{\overset{\rightarrow}{\mu}}_{2,1,3} - {\sqrt{\frac{3}{5}}c_{x}c_{z}{\overset{\rightarrow}{\mu}}_{0,0,1}} + {\sqrt{\frac{9}{5}}c_{x}{\overset{\rightarrow}{\mu}}_{1,0,2}} - {\sqrt{\frac{9}{5}}c_{z}{\overset{\rightarrow}{\mu}}_{1,1,2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,2,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\mu}}_{2,2,2}^{\prime} = 0}$${\overset{\rightarrow}{\mu}}_{2,2,3}^{\prime} = {{\overset{\rightarrow}{\mu}}_{2,2,3}^{\prime} + {\sqrt{\frac{3}{20}}\left( {c_{x}^{2} - c_{y}^{2}} \right){\overset{\rightarrow}{\mu}}_{0,0,1}} - {\sqrt{\frac{9}{5}}c_{y}{\overset{\rightarrow}{\mu}}_{1,{- 1},2}} + {\sqrt{\frac{9}{5}}c_{x}{\overset{\rightarrow}{\mu}}_{1,1,2}}}$

APPENDIX 5

First Far to Near Coefficient Transformation

The coefficients transforms simplify, and are precomputed and storedbefore running the simulation. This optimization is worth it since thereis at most 189 far interacting cells.

${\overset{\rightarrow}{\kappa}}_{0,0,0}^{\prime} = {{\frac{1}{c}{\overset{\rightarrow}{\mu}}_{0,0,1}} + {\frac{1}{{c}^{3}}{\overset{\rightarrow}{\mu}}_{0,0,3}} - {\sqrt{3}\frac{c_{y}}{{c}^{3}}{\overset{\rightarrow}{\mu}}_{1,{- 1},2}} + {\sqrt{3}\frac{c_{z}}{{c}^{3}}{\overset{\rightarrow}{\mu}}_{1,0,2}} - {\sqrt{3}\frac{c_{x}}{{c}^{3}}{\overset{\rightarrow}{\mu}}_{1,1,2}} - {\sqrt{\frac{5}{4}}\frac{{c}^{2} - {3\; c_{z}^{2}}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{2,0,3}} + {\sqrt{15}\frac{c_{x}c_{y}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{2,{- 2},3}} - {\sqrt{15}\frac{c_{y}c_{z}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{2,{- 1},3}} - {\sqrt{15}\frac{c_{x}c_{z}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{2,1,3}} + {\sqrt{\frac{15}{4}}\frac{c_{x}^{2} - c_{y}^{2}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{2,2,3}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{0,0,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{0,0,2}^{\prime} = {\frac{1}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{0,0,3}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{0,0,3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,{- 1},0}^{\prime} = 0}$${\overset{\rightarrow}{\kappa}}_{1,{- 1},1}^{\prime} = {{\sqrt{\frac{1}{3}}\frac{c_{y}}{{c}^{3}}{\overset{\rightarrow}{\mu}}_{0,0,1}} + {\sqrt{3}\frac{c_{y}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{0,0,3}} + {\frac{{c}^{2} - {3\; c_{y}^{2}}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,{- 1},2}} + {3\frac{c_{y}c_{z}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,0,2}} - {3\frac{c_{x}c_{y}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,1,2}} - {\sqrt{5}\frac{c_{x}\left( {{c}^{2} - {5\; c_{y}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,{- 2},3}} + {\sqrt{5}\frac{c_{z}\left( {{c}^{2} - {5\; c_{y}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,{- 1},3}} - {\sqrt{\frac{15}{4}}\frac{c_{y}\left( {{c}^{2} - {5\; c_{z}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,0,3}} - {\sqrt{125}\frac{c_{x}c_{y}c_{z}}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,1,3}} + {\sqrt{\frac{5}{4}}\frac{c_{y}\left( {{2{c}^{2}} + {5\left( {c_{x}^{2} - c_{y}^{2}} \right)}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,2,3}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,{- 1},2}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime} = {\sqrt{3}\frac{c_{y}}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{0,0,3}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,0,0}^{\prime} = 0}$${\overset{\rightarrow}{\kappa}}_{1,0,1}^{\prime} = {{{- \sqrt{\frac{1}{3}}}\frac{c_{z}}{{c}^{3}}{\overset{\rightarrow}{\mu}}_{0,0,1}} - {\sqrt{3}\frac{c_{z}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{0,0,3}} + {3\frac{c_{y}c_{z}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,{- 1},2}} + {\frac{{c}^{2} - {3\; c_{z}^{2}}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,0,2}} + {3\frac{c_{x}c_{z}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,1,2}} - {\sqrt{125}\frac{c_{x}c_{y}c_{z}}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,{- 2},3}} - {\sqrt{5}\frac{c_{y}\left( {{c}^{2} - {5\; c_{z}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,{- 1},3}} + {\sqrt{\frac{15}{4}}\frac{c_{z}\left( {{3{c}^{2}} - {5\; c_{z}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,0,3}} - {\sqrt{5}\frac{c_{x}\left( {{c}^{2} - {5\; c_{z}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,1,3}} - {\sqrt{\frac{125}{4}}\frac{c_{z}\left( {c_{x}^{2} - c_{y}^{2}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,2,3}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,0,2}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime} = {{- \sqrt{3}}\frac{c_{z}}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{0,0,3}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,1,0}^{\prime} = 0}$${\overset{\rightarrow}{\kappa}}_{1,1,1}^{\prime} = {{{- \sqrt{\frac{1}{3}}}\frac{c_{x}}{{c}^{3}}{\overset{\rightarrow}{\mu}}_{0,0,1}} + {\sqrt{3}\frac{c_{x}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{0,0,3}} - {3\frac{c_{x}c_{z}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,{- 1},2}} + {3\frac{c_{x}c_{z}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,0,2}} + {\frac{{c}^{2} - {3\; c_{x}^{2}}}{{c}^{5}}{\overset{\rightarrow}{\mu}}_{1,1,2}} - {\sqrt{5}\frac{c_{y}\left( {{c}^{2} - {5\; c_{x}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,{- 2},3}} - {\sqrt{125}\frac{c_{x}c_{y}c_{z}}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,{- 1},3}} - {\sqrt{\frac{15}{4}}\frac{c_{x}\left( {{c}^{2} - {5\; c_{z}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,0,3}} + {\sqrt{5}\frac{c_{z}\left( {{c}^{2} - {5\; c_{x}^{2}}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,1,3}} + {\sqrt{\frac{5}{4}}\frac{c_{x}\left( {{{- 2}{c}^{2}} + {5\left( {c_{x}^{2} - c_{y}^{2}} \right)}} \right)}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{2,2,3}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,1,2}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime} = {\sqrt{3}\frac{c_{x}}{{c}^{7}}{\overset{\rightarrow}{\mu}}_{0,0,3}}}$

APPENDIX 6

First Near Coefficient Translation

The coefficients simplify, and the first few are:

${\overset{\rightarrow}{\kappa}}_{0,0,0}^{\prime} = {{\overset{\rightarrow}{\kappa}}_{0,0,0}^{\prime} + {{c}^{2}{\overset{\rightarrow}{\kappa}}_{0,0,2}^{\prime}} - {\sqrt{3}c_{y}{\overset{\rightarrow}{\kappa}}_{1,{- 1},1}^{\prime}} - {\sqrt{3}c_{y}{c}^{2}{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime}} + {\sqrt{3}c_{z}{\overset{\rightarrow}{\kappa}}_{1,0,1}^{\prime}} + {\sqrt{3}c_{z}{c}^{2}{\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime}} - {\sqrt{3}c_{x}{\overset{\rightarrow}{\kappa}}_{1,1,1}^{\prime}} - {\sqrt{3}c_{x}{c}^{2}{\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime}} + {\sqrt{15}c_{x}c_{y}{\overset{\rightarrow}{\kappa}}_{2,{- 2},2}^{\prime}} - {\sqrt{15}c_{y}c_{z}{\overset{\rightarrow}{\kappa}}_{2,{- 1},2}^{\prime}} - {\sqrt{\frac{5}{4}}\left( {{c}^{2} - {3\; c_{z}^{2}}} \right){\overset{\rightarrow}{\kappa}}_{2,0,2}^{\prime}} - {\sqrt{15}c_{x}c_{z}{\overset{\rightarrow}{\kappa}}_{2,1,2}^{\prime}} + {\sqrt{\frac{15}{4}}\left( {c_{x}^{2} - c_{y}^{2}} \right){\overset{\rightarrow}{\kappa}}_{2,2,2}^{\prime}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{0,0,1}^{\prime} = 0}$${\overset{\rightarrow}{\kappa}}_{0,0,2}^{\prime} = {{\overset{\rightarrow}{\kappa}}_{0,0,2} - {\sqrt{\frac{25}{3}}c_{y}{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}} + {\sqrt{\frac{25}{3}}c_{z}{\overset{\rightarrow}{\kappa}}_{1,0,3}} - {\sqrt{\frac{25}{3}}c_{x}{\overset{\rightarrow}{\kappa}}_{1,1,3}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{0,0,3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,{- 1},0}^{\prime} = 0}$${\overset{\rightarrow}{\kappa}}_{1,{- 1},1}^{\prime} = {{{- \sqrt{\frac{4}{3}}}c_{y}{\overset{\rightarrow}{\kappa}}_{0,0,2}^{\prime}} + {\overset{\rightarrow}{\kappa}}_{1,{- 1},1}^{\prime} + {\left( {{c}^{2} - {5\; c_{y}^{2}}} \right){\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime}} - {2\; c_{y}c_{z}{\overset{\rightarrow}{\kappa}}_{1,0,3}} + {2\; c_{x}c_{y}{\overset{\rightarrow}{\kappa}}_{1,1,3}} - {\sqrt{5}c_{x}{\overset{\rightarrow}{\kappa}}_{2,{- 2},2}^{\prime}} + {\sqrt{5}c_{z}{\overset{\rightarrow}{\kappa}}_{2,{- 1},2}^{\prime}} + {\sqrt{\frac{5}{3}}c_{y}{\overset{\rightarrow}{\kappa}}_{2,0,2}^{\prime}} + {\sqrt{5}c_{y}{\overset{\rightarrow}{\kappa}}_{2,2,2}^{\prime}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,{- 1},2}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime} = {\overset{\rightarrow}{\kappa}}_{1,{- 1},3}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,0,0}^{\prime} = 0}$${\overset{\rightarrow}{\kappa}}_{1,0,1}^{\prime} = {{\sqrt{\frac{4}{3}}c_{z}{\overset{\rightarrow}{\kappa}}_{0,0,2}^{\prime}} - {2\; c_{y}c_{z}{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime}} + {\overset{\rightarrow}{\kappa}}_{1,0,1}^{\prime} + {\left( {{c}^{2} + {2\; c_{z}^{2}}} \right){\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime}} - {2\; c_{x}c_{z}{\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime}} - {\sqrt{5}c_{y}{\overset{\rightarrow}{\kappa}}_{2,{- 1},2}^{\prime}} + {\sqrt{\frac{20}{3}}c_{z}{\overset{\rightarrow}{\kappa}}_{2,0,2}^{\prime}} - {\sqrt{5}c_{x}{\overset{\rightarrow}{\kappa}}_{2,1,2}^{\prime}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,0,2}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime} = {\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,1,0}^{\prime} = 0}$${\overset{\rightarrow}{\kappa}}_{1,1,1}^{\prime} = {{{- \sqrt{\frac{4}{3}}}c_{x}{\overset{\rightarrow}{\kappa}}_{0,0,2}^{\prime}} + {2\; c_{x}c_{y}{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime}} - {2\; c_{x}c_{z}{\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime}} + {\overset{\rightarrow}{\kappa}}_{1,1,1}^{\prime} + {\left( {{c}^{2} + {2\; c_{x}^{2}}} \right){\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime}} - {\sqrt{5}c_{y}{\overset{\rightarrow}{\kappa}}_{2,{- 2},2}^{\prime}} + {\sqrt{\frac{5}{3}}c_{x}{\overset{\rightarrow}{\kappa}}_{2,0,2}^{\prime}} + {\sqrt{5}c_{z}{\overset{\rightarrow}{\kappa}}_{2,1,2}^{\prime}} - {\sqrt{5}c_{x}{\overset{\rightarrow}{\kappa}}_{2,2,2}^{\prime}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,1,2}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime} = {\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,{- 2},0}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,{- 2},1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,{- 2},2}^{\prime} = {{{- \sqrt{\frac{4}{5}}}c_{x}{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime}} - {\sqrt{\frac{4}{5}}c_{y}{\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime}} + {\overset{\rightarrow}{\kappa}}_{2,{- 2},2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,{- 2},3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,{- 1},0}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,{- 1},1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,{- 1},2}^{\prime} = {{{- \sqrt{\frac{4}{5}}}c_{z}{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime}} - {\sqrt{\frac{4}{5}}c_{y}{\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime}} + {\overset{\rightarrow}{\kappa}}_{2,{- 1},2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,{- 1},3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,0,0}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,0,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,0,2}^{\prime} = {{\sqrt{\frac{4}{15}}c_{y}{\overset{\rightarrow}{\kappa}}_{1,{- 1},3}^{\prime}} - {\sqrt{\frac{16}{15}}c_{z}{\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime}} + {\sqrt{\frac{4}{15}}c_{x}{\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime}} + {\overset{\rightarrow}{\kappa}}_{2,0,2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,0,3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,1,0}^{\prime} = 0}$$\mspace{79mu} {{\overset{\rightarrow}{\kappa}}_{2,1,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,1,2}^{\prime} = {{{- \sqrt{\frac{4}{5}}}c_{x}{\overset{\rightarrow}{\kappa}}_{1,0,3}^{\prime}} + {\sqrt{\frac{4}{5}}c_{z}{\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime}} + {\overset{\rightarrow}{\kappa}}_{2,1,2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,1.3}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,2,0}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,2,1}^{\prime} = 0}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,2,2}^{\prime} = {{\sqrt{\frac{4}{5}}c_{y}{\overset{\rightarrow}{\kappa}}_{1,{- {,3}}}^{\prime}} - {\sqrt{\frac{4}{5}}c_{x}{\overset{\rightarrow}{\kappa}}_{1,1,3}^{\prime}} + {\overset{\rightarrow}{\kappa}}_{2,2,2}}}$$\mspace{20mu} {{\overset{\rightarrow}{\kappa}}_{2,2,3}^{\prime} = 0}$

APPENDIX 7

Buoyancy

The vorticity induced by buoyancy can be computed directly by samplingeverywhere in R³ the buoyancy term with vorticles. But this laboriousintegral can be avoided, since buoyancy vorticity is concentrated nearplaces of varying density and where the density gradient and buoyantacceleration are perpendicular. First, a set of vorticles can beassociated with a single particle j, which can be generalized tomultiple density particles. Let us define a local

coordinate system centered at x_(j) such that {right arrow over (z)} isaligned with

$\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - {\frac{d\; \overset{\rightarrow}{u}}{dt}.}$

For a density particle or radius r_(j), it can be expected that thevorticles induced by buoyancy are concentrated around a ring {x_(α),αϵ2π} of diameter r_(j), perpendicular to vector unit {right arrow over(z)}, with an axis {right arrow over (n)}_(α). In local coordinates:

$\begin{matrix}{{x_{\alpha} = {\frac{r_{j}}{2}\left( {{\cos (\alpha)},{\sin (\alpha)},0} \right)}}{{\overset{\rightarrow}{n}}_{\alpha} = {\left( {{- {\sin (\alpha)}},{\cos (\alpha)},0} \right).}}} & (52)\end{matrix}$

For a canonical density particle, the following density field may beused:

$\begin{matrix}{{\hat{\rho}}_{j} = {\exp\left( \left( {1 + {\kappa_{1}\frac{x \cdot x}{2r_{j}^{2}}}} \right)^{- \kappa_{2}} \right)}} & (53)\end{matrix}$

The parameters κ₀, κ₁ and κ₂ may be fitted using the following equation:

∇ log({circumflex over (ρ)}_(j))×{right arrow over (z)}=r _(j)²κ₀∇×(∫_(α)ξ_(i)(p−x _(α))³ {right arrow over (n)} _(α)×(p−x _(α))dα).  (54)

The following parameters may be obtained:

κ₀=0.493483 κ₁=0.572636 κ₂=3.423340  (55)

FIG. 26A shows two almost exactly overlapping curves of the largestcomponent of the left (red) and right (blue) side of Eq. (54) accordingto some embodiments. FIG. 26B shows the density field ρ_(j) for r_(j)=1and m_(j)=1. The density fits the closed form integral with a thirddecimal accuracy on all 3 vector components.

To generalize to multiple particles, {circumflex over (ρ)}_(j) can beremapped to values in (0, m_(j)) so that densities can be added, and thefield of a single particle density

$\rho_{j} = {m_{j}\frac{{\hat{\rho}}_{j} - 1}{e - 1}}$

can be obtained. Since

${{\frac{\nabla\rho}{\rho} \times \overset{\rightarrow}{z}} = {\sum{\frac{m_{j}}{e - 1}\frac{{\hat{\rho}}_{j}}{\rho}{\nabla{\log \left( {\hat{\rho}}_{j} \right)}} \times \overset{\rightarrow}{z}}}},$

a set of weights can be obtained that modulate the rotations of thevorticles induced by the density particle j, that accounts foracceleration and multiple particles:

$\begin{matrix}{{\overset{\rightarrow}{w}}_{\alpha} = {r_{j}^{2}\kappa_{0}\frac{\pi}{n}{{\overset{\rightarrow}{g} + \overset{\rightarrow}{e} - \frac{d\overset{\rightarrow}{u}}{dt}}}\frac{m_{j}}{e - 1}\frac{{\hat{\rho}}_{j}\left( x_{\alpha} \right)}{\rho \left( x_{\alpha} \right)}{\overset{\rightarrow}{n}}_{\alpha}}} & (56)\end{matrix}$

APPENDIX 8

Boundary Vorticity

The vorticity near the boundary δΩ of a moving object is given by∇log(ρ)×({right arrow over (F)}−{right arrow over (a)})+∇×{right arrowover (F)}, where {right arrow over (a)} is the acceleration near theboundary. To compute this term, one may consider a coordinate systemwhere the half space z<0 is inside the solid object, z=0 in on theobject's surface, and z>0 is outside. Using the Heavyside step functionH , the nearby density is formally defined as ρ=H(−z)ρ₀+H(z)ρ₁, theexternal forces term is defined as {right arrow over (F)}=H(−z){rightarrow over (f)}+H(z)({right arrow over (g)}+{right arrow over (e)}), andthe acceleration is

$\overset{\rightarrow}{a} = {{{H\left( {- z} \right)}\overset{\rightarrow}{f}} + {{H(z)}{\frac{d\overset{\rightarrow}{u}}{dt}.}}}$

Since vorticity is solved for, {right arrow over (F)} can be replacedwith any {right arrow over (F)}+{right arrow over (G)} to satisfy theboundary condition as long as ∇×{right arrow over (G)}=0 outside ofsolid objects. First let us consider the following {right arrow over(G)}:

$\begin{matrix}{\overset{\rightarrow}{G} = {{{H(z)}\left( {\frac{d\overset{\rightarrow}{u}}{dt} - \overset{\rightarrow}{g} - {\frac{1}{2}\overset{\rightarrow}{e}}} \right)} - {\frac{1}{2}{\overset{\rightarrow}{e}.}}}} & (57)\end{matrix}$

Although

$\frac{d\overset{\rightarrow}{u}}{dt}$

may have curl, the object surface may be discretized in regions ofconstant

$\frac{d\overset{\rightarrow}{u}}{dt}.$

To prevent the space inside objects from participating to the non-rigidfluid motion, one may set ρ₀→∞, thus enforcing {right arrow over (f)}inside objects. One may take the limit of the integral of the collisionterm near the surface, using the following limit to the Heavyside stepfunction

${{H(z)} = {{\lim\limits_{t->\infty}\frac{1}{2}} + {\frac{1}{\pi}{\tan^{- 1}({zt})}}}},$

and obtain the surface vorticity:

$\begin{matrix}{{{\lim\limits_{\rho_{0}->\infty}{\lim\limits_{t->\infty}{\int_{- h}^{h}{{\nabla{\log (\rho)}} \times \left( {\overset{\rightarrow}{F} + \overset{\rightarrow}{G} - \frac{d\overset{\rightarrow}{u}}{dt}} \right)}}}} + {\nabla{\times \left( {\overset{\rightarrow}{F} + \overset{\rightarrow}{G}} \right){dz}}}} = {\left( {\overset{\rightarrow}{f} - \overset{\rightarrow}{e} - \overset{\rightarrow}{a}} \right) \times \overset{\rightarrow}{n}}} & (58)\end{matrix}$

APPENDIX 9

Source Fields

The source integral of Eq. (26) is given by van Oosterom [2011]. Beloware the formula of the source panel field S, and its derivative ∇S, fora triangle defined by {p₀, p₁, p₂}, using the symbols of Eq. (60):

$\begin{matrix}{{{S(p)} = {\frac{1}{\overset{\rightarrow}{n}}\left( {{\frac{\overset{\rightarrow}{n} \cdot \left( {{\overset{\rightarrow}{d}}_{1} \times {\overset{\rightarrow}{d}}_{2}} \right)}{e_{0}}{\log\left( \frac{f_{0}}{g_{0}} \right)}} + {\frac{\overset{\rightarrow}{n} \cdot \left( {{\overset{\rightarrow}{d}}_{2} \times {\overset{\rightarrow}{d}}_{0}} \right)}{e_{1}}{\log\left( \frac{f_{1}}{g_{1}} \right)}} + {\frac{\overset{\rightarrow}{n} \cdot \left( {{\overset{\rightarrow}{d}}_{0} \times {\overset{\rightarrow}{d}}_{1}} \right)}{e_{2}}{\log\left( \frac{f_{2}}{g_{2}} \right)}} - {2k_{1}a\; \tan \; 2\left( {k_{1},k_{0}} \right)}} \right)}}{{\nabla{S(p)}} = {\frac{1}{\overset{\rightarrow}{n}}\left( {{\left( {\overset{\rightarrow}{n} \cdot \left( {{\overset{\rightarrow}{d}}_{1} \times {\overset{\rightarrow}{d}}_{2}} \right)} \right){\overset{\rightarrow}{h}}_{0}} + {\frac{\left( {{\overset{\rightarrow}{d}}_{1} - {\overset{\rightarrow}{d}}_{2}} \right) \times \overset{\rightarrow}{n}}{e_{0}}{\log\left( \frac{f\; 0}{g\; 0} \right)}} + {\left( {\overset{\rightarrow}{n} \cdot \left( {{\overset{\rightarrow}{d}}_{2} \times {\overset{\rightarrow}{d}}_{0}} \right)} \right){\overset{\rightarrow}{h}}_{1}} + {\frac{\left( {{\overset{\rightarrow}{d}}_{2} - {\overset{\rightarrow}{d}}_{0}} \right) \times \overset{\rightarrow}{n}}{e_{1}}{\log\left( \frac{f\; 1}{g\; 1} \right)}} + {\left( {\overset{\rightarrow}{n} \cdot \left( {{\overset{\rightarrow}{d}}_{0} \times {\overset{\rightarrow}{d}}_{1}} \right)} \right){\overset{\rightarrow}{h}}_{2}} + {\frac{\left( {{\overset{\rightarrow}{d}}_{0} - {\overset{\rightarrow}{d}}_{1}} \right) \times \overset{\rightarrow}{n}}{e_{2}}{\log\left( \frac{f\; 2}{g\; 2} \right)}} + {{D(p)}\overset{\rightarrow}{n}} + {k_{1}{\nabla{D(p)}}}} \right)}}} & (59)\end{matrix}$

The symbols above are defined as:

$\begin{matrix}{{{\overset{\rightarrow}{n} = {\left( {p_{1} - p_{0}} \right) \times \left( {p_{2} - p_{0}} \right)}}{{\overset{\rightarrow}{d}}_{i} = {p_{i} - p}}{l_{i} = {{\overset{\rightarrow}{d}}_{i}}}{k_{0} = {{l_{0}l_{1}l_{2}} + {l_{0}{{\overset{\rightarrow}{d}}_{1} \cdot {\overset{\rightarrow}{d}}_{2}}} + {l_{1}{{\overset{\rightarrow}{d}}_{2} \cdot {\overset{\rightarrow}{d}}_{0}}} + {l_{2}{{\overset{\rightarrow}{d}}_{0} \cdot {\overset{\rightarrow}{d}}_{1}}}}}{k_{1} = {\left( {{\overset{\rightarrow}{d}}_{2} \times {\overset{\rightarrow}{d}}_{1}} \right) \cdot {\overset{\rightarrow}{d}}_{0}}}{e_{i} = {{p_{j} - p_{k}}}}{f_{i} = {l_{j} + l_{k} + e_{i}}}{g_{i} = {l_{j} + l_{k} - e_{i}}}{{\overset{\rightarrow}{h}}_{i} = {\frac{2}{f_{i}g_{i}}\left( {\frac{{\overset{\rightarrow}{d}}_{j}}{l_{j}} + \frac{{\overset{\rightarrow}{d}}_{k}}{l_{k}}} \right)}}j = {\left( {i + 1} \right)\mspace{14mu} {\% 3}}}{k = {\left( {i + 2} \right)\mspace{14mu} {\% 3}}}} & (60)\end{matrix}$

The function ∇S(p) is numerically unstable when p is aligned with theedges. This corresponds to a situation where some triangles in thedecomposition of Carley [2012] are very thin. FIGS. 27A-27D illustratethat the integral over any triangle can be decomposed into 6 right-angletriangle integrals. FIG. 27A shows that point p is the projection of theevaluation point p onto the plane of the triangle. FIG. 27B shows thatthe triangle is split into three signed triangles with a corner at p.FIG. 27C shows that each triangle is split into two signed right angletriangles. FIG. 27D shows that the integral of a piece is evaluated withp as the origin and an axis aligned triangle.

Using this decomposition, the thin triangles may be eliminated. Thesource field induced at a point p located on the z axis by a trianglewith a right angle at p₀ lying in the plane z=0 with corner p₂ locatedat the origin is given by:

∇{tilde over (S)}(p)_(x)=arctanh(k ₁ /k ₂)+k ₁ k ₇

∇{tilde over (S)}(p)_(y) =k ₆−log(1+k ₅)−k ₇

∇{tilde over (S)}(p)_(z)=atan2(k ₁(k ₂ −k ₀), k ₀ k ₁ ² +k ₂)  (61)

The symbols above are defined as:

$\begin{matrix}{{k_{0} = {{p}/{p_{0}}}}{k_{1} = {{{p_{1} - p_{0}}}/{p_{0}}}}{k_{2} = \sqrt{1 + k_{1}^{2} + k_{0}^{2}}}{k_{3} = {{\left( {k_{1}^{2} - 1} \right)\left( {k_{0}^{2} - 1} \right)} - 2}}{k_{4} = \sqrt{1 + k_{1}^{2}}}{k_{5} = \sqrt{1 + k_{0}^{2}}}{k_{6} = {\log \left( k_{0} \right)}}{k_{7} = {\left( {k_{6} - {\log \left( {k_{2} + k_{4}} \right)}} \right)/k_{4}}}} & (62)\end{matrix}$

APPENDIX 10

Doublet Fields

The doublet integral of Eq. (26) is given by Oosterom and Strackee[1983]. Below are the formula of the doublet panel field D, and itsderivative ∇D, for a triangle defined by {p₀, p₁,p₂}, using the symbolsof Eq. (64):

$\begin{matrix}{{{D(p)} = {{- 2}a\; \tan \; 2\left( {k_{1},k_{0}} \right)}}{{\nabla{D(p)}} = {{- 2}\frac{{k_{0}\overset{\rightarrow}{n}} - {k_{1}{\overset{\rightarrow}{k}}_{2}}}{k_{0}^{2} + k_{1}^{2}}}}} & (63)\end{matrix}$

The symbols above are defined as:

$\begin{matrix}{\mspace{79mu} {{\overset{\rightarrow}{n} = {\left( {p_{1} - p_{0}} \right) \times \left( {p_{2} - p_{0}} \right)}}\mspace{20mu} {{\overset{\rightarrow}{d}}_{i} = {p_{i} - p}}\mspace{20mu} {l_{i} = {{\overset{\rightarrow}{d}}_{i}}}\mspace{20mu} {k_{0} = {{l_{0}l_{1}l_{2}} + {l_{0}{{\overset{\rightarrow}{d}}_{1} \cdot {\overset{\rightarrow}{d}}_{2}}} + {l_{1}{{\overset{\rightarrow}{d}}_{2} \cdot {\overset{\rightarrow}{d}}_{0}}} + {l_{2}{{\overset{\rightarrow}{d}}_{0} \cdot {\overset{\rightarrow}{d}}_{1}}}}}\mspace{20mu} {k_{1} = {\left( {{\overset{\rightarrow}{d}}_{2} \times {\overset{\rightarrow}{d}}_{1}} \right) \cdot {\overset{\rightarrow}{d}}_{0}}}{{\overset{\rightarrow}{k}}_{2} = {{{- \left( {\frac{{l_{1}l_{2}} + {{\overset{\rightarrow}{d}}_{1} \cdot {\overset{\rightarrow}{d}}_{2}}}{l_{0}} + l_{1} + l_{2}} \right)}{\overset{\rightarrow}{d}}_{0}} - {\left( {\frac{{l_{0}l_{2}} + {{\overset{\rightarrow}{d}}_{2} \cdot {\overset{\rightarrow}{d}}_{0}}}{l_{1}} + l_{2} + l_{0}} \right){\overset{\rightarrow}{d}}_{1}} - {\left( {\frac{{l_{0}l_{1}} + {{\overset{\rightarrow}{d}}_{0} \cdot {\overset{\rightarrow}{d}}_{1}}}{l_{2}} + l_{0} + l_{1}} \right){\overset{\rightarrow}{d}}_{2}}}}}} & (64)\end{matrix}$

What is claimed is:
 1. A method of computer graphic image processing,wherein a computer system generates rendered image frames ofincompressible gas corresponding to a camera view of a three-dimensionalvirtual space, the method comprising: generating a plurality ofvorticles, each respective vorticle centered at a respective position inthe virtual space and having a respective spatial size; computingadvected field of the plurality of vorticles by solving a vorticityequation using a fast multipole method (FMM); dividing the virtual spaceinto a plurality of concentric regions surrounding a position of acamera; defining a respective density voxel grid for each respectiveregion of the plurality of concentric regions, each respective densityvoxel grid having a respective first voxel size that increases withincreasing distance from the position of the camera; defining arespective velocity voxel grid for each respective region of theplurality of concentric regions, each respective velocity voxel gridhaving a respective second voxel size that increases with increasingdistance from the position of the camera; adding a new densityrepresenting a smoke emitter, wherein the new density has a firstdistribution in a first region of the virtual space, wherein the firstregion overlaps with one or more density voxels in one or more regionsof the plurality of concentric regions; activating and adding density toeach of the one or more density voxels of the one or more regions thatoverlap with the first region; activating and computing velocity fieldat each velocity voxel of the one or more regions that overlaps with theactive density voxels; moving the new density based on the velocityfield of each active velocity voxel to obtain a second distribution ofthe new density in a second region of the virtual space for a nextframe; activating and adding density to each density voxel that overlapswith the second region; and storing the density of each active densityvoxel and the velocity field of each active velocity voxel in acomputer-readable medium as the next frame for display at a screen. 2.The method of claim 1, wherein the respective second voxel size of eachrespective velocity voxel grid is greater than the respective firstvoxel size of the corresponding density voxel grid.
 3. The method ofclaim 1, wherein the plurality of concentric regions is defined by aplurality of concentric spheres, each respective region located betweentwo consecutive spheres.
 4. The method of claim 1, wherein computing theadvected field of the plurality of vorticles using the FMM comprises:building an octree of vorticles including a plurality octree levels,each respective octree level having a cell size proportional to 2^(−L),L being a level number of the respective octree level; storing eachrespective vorticle of the plurality of vorticles in a respective cellof a respective octree level having a cell size that is proportional tothe spatial size of the respective vorticle; allocating near cells as afirst cell containing the position of the camera and cells that aredirect neighbors of the first cell, and far interacting cells as allother cells; and computing the advected field using the FMM by summingthe vorticles stored in each respective far cell.
 5. The method of claim4, wherein each respective vorticle is stored in the respective cellhaving the cell size of$2^{\frac{\log \mspace{11mu} r_{i}}{\log \mspace{11mu} 2}},$ r_(i)being the spatial size of the respective vorticle.
 6. The method ofclaim 4, further comprising: computing velocity field of the pluralityof vorticles; and moving the plurality of vorticles to a next framebased on the velocity field using the octree.
 7. The method of claim 4,wherein the respective spatial size of each respective vorticleincreases with increasing distance between the respective position ofthe respective vorticle and the position of the camera.
 8. The method ofclaim 1, further comprising: inputting an object defined by a boundarysurface; dividing the boundary surface into a plurality of panels, eachrespective panel centered at a respective position and having arespective panel size; computing advected field at each respective panelusing the FMM; computing a source panel coefficient and a doublet panelcoefficient for each respective panel; and computing harmonic field bysolving the vorticity equation using the FMM based on the source panelcoefficients and the doublet panel coefficients of the plurality ofpanels.
 9. The method of claim 8, wherein computing the harmonic fieldusing the FMM comprises: building an octree of panels including aplurality octree levels, each respective octree level having a cell sizeproportional to 2^(−L), L being a level number of the respective octreelevel; storing each respective panel of the plurality of panels in arespective cell of a respective octree level having a cell size that isproportional to the spatial size of the respective panel; allocatingnear cells as a first cell containing the position of the camera andcells that are direct neighbors of the first cell, and far interactingcells as all other cells; and computing the harmonic field using the FMMby summing the panels stored in each respective far cell.
 10. The methodof claim 1, further comprising: stretching one or more vorticles of theplurality of vorticles by a stretching factor s; and unstretching theone or more vorticles while preserving an energy and enstrophy of eachof the one or more vorticles.
 11. A computer product comprising anon-transitory computer readable medium storing a plurality ofinstructions that when executed control a computer system to generaterendered image frames of incompressible gas corresponding to a cameraview of a three-dimensional virtual space, wherein the plurality ofinstructions comprises: generating a plurality of vorticles, eachrespective vorticle centered at a respective position in the virtualspace and having a respective spatial size; computing advected field ofthe plurality of vorticles by solving a vorticity equation using a fastmultipole method (FMM); dividing the virtual space into a plurality ofconcentric regions surrounding a position of a camera; defining arespective density voxel grid for each respective region of theplurality of concentric regions, each respective density voxel gridhaving a respective first voxel size that increases with increasingdistance from the position of the camera; defining a respective velocityvoxel grid for each respective region of the plurality of concentricregions, each respective velocity voxel grid having a respective secondvoxel size that increases with increasing distance from the position ofthe camera; adding a new density representing a smoke emitter, whereinthe new density has a first distribution in a first region of thevirtual space, wherein the first region overlaps with one or moreregions of the plurality of concentric regions; activating and addingdensity to each density voxel of the one or more regions that overlapswith the first region; activating and computing velocity field at eachvelocity voxel of the one or more regions that overlaps with the activedensity voxels; moving the new density based on the velocity field ofeach active velocity voxel to obtain a second distribution of the newdensity in a second region of the virtual space for a next frame;activating and adding density to each density voxel that overlaps withthe second region; and storing the density of each active density voxeland the velocity field of each active velocity voxel in acomputer-readable medium as the next frame for display at a screen. 12.The computer product of claim 11, wherein the respective second voxelsize of each respective velocity voxel grid is greater than therespective first voxel size of the corresponding density voxel grid. 13.The computer product of claim 11, wherein the plurality of concentricregions is defined by a plurality of concentric spheres, each respectiveregion located between two consecutive spheres.
 14. The computer productof claim 11, wherein the instructions for computing the advected fieldof the plurality of vorticles using the FMM comprise: building an octreeof vorticles including a plurality octree levels, each respective octreelevel having a cell size proportional to 2^(−L), L being a level numberof the respective octree level; storing each respective vorticle of theplurality of vorticles in a respective cell of a respective octree levelhaving a cell size that is proportional to the spatial size of therespective vorticle; allocating near cells as a first cell containingthe position of the camera and cells that are direct neighbors of thefirst cell, and far interacting cells as all other cells; and computingthe advected field using the FMM by summing the vorticles stored in eachrespective far cell.
 15. The computer product of claim 14, wherein eachrespective vorticle is stored in the respective cell having the cellsize of $2^{\frac{\log \mspace{11mu} r_{i}}{\log \mspace{11mu} 2}},$r_(i) being the spatial size of the respective vorticle.
 16. Thecomputer product of claim 14, wherein the plurality of instructionsfurther comprises: computing velocity field of the plurality ofvorticles; and moving the plurality of vorticles to a next frame basedon the velocity field using the octree.
 17. The computer product ofclaim 14, wherein the respective spatial size of each respectivevorticle increases with increasing distance between the respectiveposition of the respective vorticle and the position of the camera. 18.The computer product of claim 11, wherein the plurality of instructionsfurther comprises: inputting an object defined by a boundary surface;dividing the boundary surface into a plurality of panels, eachrespective panel centered at a respective position and having arespective panel size; computing advected field at each respective panelusing the FMM; computing a source panel coefficient and a doublet panelcoefficient for each respective panel; and computing harmonic field bysolving the vorticity equation using the FMM based on the source panelcoefficients and the doublet panel coefficients of the plurality ofpanels.
 19. The computer product of claim 18, wherein the instructionsfor computing the harmonic field using the FMM comprise: building anoctree of panels including a plurality octree levels, each respectiveoctree level having a cell size proportional to 2^(−L), L being a levelnumber of the respective octree level; storing each respective panel ofthe plurality of panels in a respective cell of a respective octreelevel having a cell size that is proportional to the spatial size of therespective panel; allocating near cells as a first cell containing theposition of the camera and cells that are direct neighbors of the firstcell, and far interacting cells as all other cells; and computing theharmonic field using the FMM by summing the panels stored in eachrespective far cell.
 20. The computer product of claim 11, wherein theplurality of instructions further comprises: stretching one or morevorticles of the plurality of vorticles by a stretching factor s; andunstretching the one or more vorticles while preserving an energy andenstrophy of each of the one or more vorticles.